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Chapter  I 


IMTRODUC  TION 

The  low-density,  high-speed  flow  of  a  viscous  fluid  around  a  blunt 
body  presents  a  wide  variety  of  problems  that  have  attracted  considerable 
attention.  One  such  area  of  interest  is  concerned  with  the  effects  of  a 
rarefaction  of  the  fluid  at  high  altitudes.  This  rarefaction  can  lead 
to  a  wide  range  of  flow  conditions  --  from  continuum  flow  at  low  altitudes 
to  free -molecular  flow  at  very  high  altitudes.  The  transition  region 
between  these  two  limits  has  been  subdivided  into  various  flow  regimes 
depending  on  the  degree  of  rarefaction  of  the  fluid.  These  regimes  have 
been  described  by  Hayes  and  Probstein  [l8]*,  Frobstein  [27],  and  Cheng  [2]. 
These  regimes  are  not  sharply  defined,  and  most  investigations  have  been 
concerned  with  a  penetration  of  the  transition  region  from  either  the 
continuum  or  free-molecule  end.  The  present  investigation  will  be 
concerned  with  just  such  a  penetration  from  the  continuum  end* 

In  the  basic  continuum  approach,  valid  for  large  Reynolds  numbers, 
the  flow  between  the  detached  shock  and  the  body  is  divided  into  two 
regions:  an  outer  region  in  which  the  viscous  forces  are  negligible 
extends  across  most  of  the  shock  layer;  and  a  narrow  inner  region  that 
can  be  treated  by  Prandtl's  boundary -layer  theory  contains  all  the 
effects  of  viscosity.  At  high  altitudes  where  the  Reynolds  number  is 


The  numbers  enclosed  in  brackets  refer  to  the  list  of  references  at 
the  end  of  the  paper. 
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not  large,  the  viscous  layer  thickens  and  the  classical  houndary-layer 
theory  is  no  longer  adequate.  This  flow  regime  —  in  which  the  viscous 
layer  is  too  thick  to  be  treated  by  classical  boundary- layer  theory 
but  is  still  smaller  than  the  shock-layer  thickness,  and  in  which  the 
effect  of  a  thickening  of  the  shock  is  not  appreciable  —  is  the  region 
treated  in  this  paper. 

The  Navier-Stokes  equations,  modified  by  appropriate  slip  and 
temperature-, lump  conditions  at  the  boundaries,  form  the  basis  of 
such  a  study.  The  study  of  these  equations  for  low-density  flow  has 
been  the  object  of  numerous  investigations,  and  several  methods  of 
approach  have  been  developed.  Of  primary  interest  for  the  present 
investigation  are  the  second-order  boundary- layer  theory,  the  method 
of  local  similarity  (and  a  generalization  --  the  method  of  series 
truncation),  and  the  thin-shock-layer  analysis  of  H.K.  Cheng  [1,2,5] • 

The  classical  boundary- layer  theory  of  Prandtl  yields  asymptotic 
solutions,  valid  at  large  Reynolds  number,  to  the  Navier-Stokes 
equations.  The  analysis  of  secondary  boundary- layer  effects,  which 
become  important  at  lower  Reynolds  numbers,  has  been  obtained  by 
application  of  the  method  of  matched  asymptotic  expansions  to  the 
Navier-Stokes  equations.  The  resulting  second-order  boundary -layer 
theory,  which  has  been  developed  most  systematically  by  Van  Dyke  [52], 
represents  an  extension  of  continuum  results  toward  the  regime  of 
rarefied  flow.  The  secondary  effects  which  are  analysed  in  this  manner 
are  the  effects  of  the  curvature  of  the  body,  the  slip  and  temperature - 
jump  conditions  at  the  boundary,  the  vorticity  in  the  outer  flow,  and 


the  perturbation  of  the  outer  flow  by  the  displacement  thickness  of 
the  first-order  boundary  layer.  A  discussion  of  the  development  of 
second-order  boundary- layer  theory  up  to  1962  is  contained  in 
reference  33*  In  addition,  extensive  computation  of  the  second-order 
boundary-layer  effects  have  been  carried  out  by  Davis  and  Flugge-Lotz 
[12]  and  Fannelop  and  Flugge-Lotz  [  13l j?  using  an  implicit  finite- 
difference  method. 

This  second-order  boundary-layer  theory  has  yielded  results  of 
considerable  interest,  but  some  practical  difficulties  remain  in  its 
application.  Its  use  requires  the  solution  of  a  series  of  inter¬ 
related  problems,  each  of  which  is  simpler  in  concept  than  the  basic 
problem  defined  by  the  Navier-Stokes  equations.  However,  despite  the 
apparent  simplicity  of  the  inviscid,  outer  flow  (on  which  the  boundary- 
layer  solutions  are  based),  the  available  methods  of  solution  for  the 
inviscid  flow  are  frequently  inadequate.  The  simpler,  analytical 
methods  of  solution  often  encounter  convergence  difficulties  [37),  and 
even  the  more  elaborate  numerical  procedures  may  lean  to  results  of 
insufficient  reliability  [25 ]•  Due  to  the  difficulty  of  obtaining  such 
inviscid  solutions,  many  investigators  have  used  various  approximations 
to  determine  the  second-order  outer  flow  due  to  displacement  thickness, 
and  this  affects  the  accuracy  of  their  results . 

As  an  alternative,  Davis  and  Flugge-Lotz  proposed  that  the  flow 
be  represented  by  a  simplified  form  of  the  Navier-Stokes  equations 
that  would  be  uniformly  valid  throughout  the  shock  layer.  The  equations 
would  contain  only  those  terms  that  contribute  to  the  first-  or 
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second-order  boundary- layer  equations  or  to  the  first-  or  second-order 
outer,  inviscid  flow.  Thus  the  equations  are  valid  in  the  same  flow 
regime  as  the  second-order  boundary-layer  theory,  but  they  present  a 
unified  treatment  of  the  entire  shock  layer.  The  feature  of  these 
simplified  equations  that  makes  their  solutions  more  feasible  than  the 
solution  of  the  entire  Navier-Stokes  equations  is  that  all  the 
characteristic  surfaces  of  the  simplified  equations  are  real.  This 
feature  arises  the  possibility  of  solving  the  equations  as  an  initial- 
value  problem1,  and  the  purpose  of  the  present  investigation  is  to 
examine  methods  of  analysis  that  lead  to  such  a  solution.  Two  general 
methods  of  analysis  are  considered-  First,  we  obtain  solutions  of  the 
equations  that  are  valid  in  the  stagnation  region  of  the  flow.  Second, 
we  examine  the  use  of  implicit  finite-difference  methods  to  solve  the 
equations  as  an  initial-value  problem  using  the  solutions  at  the  axis 
as  initial  data. 

The  method  of  local  similarity  has  been  widely  employed  in  the 
study  of  low-density  flow  in  the  neighborhood  of  the  axis  of  symmetry 
of  a  blunt  body.  In  this  method,  an  assumed  similitude  is  used  to 
effect  a  separation  of  variables,  and  the  equations  are  then  integrated 
along  the  stagnation  streamline.  In  addition,  as  a  numerical  convenience, 
the  basic  differential  equations  are  often  simplified  by  the  introduction 


The  problem  is  more  accurately  described  as  an  initial-boundary-value 
problem.  These  two  descriptions  are  used  somewhat  interchangeably  in 
this  paper,  bu*-  this  should  not  present  any  difficulty. 
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of  additional  approximations.  The  approximations  which  have  been  most 
commonly  used  are  a  constant -density  flow,  e.g.  Probstein  and  Kemp 
[28],  and  the  thin-shock-layer  model,  e.g.  Ho  and  ProbBtein  [19]- 
The  assumption  of  local  similitude  is  based  on  the  existence  of  a 
spherical  symmetry  in  the  geometry  of  the  flow  boundaries,  and  the 
thin-shock- layer  concept  is  frequently  used  as  Justification  for 
this  latter  assumption. 

A  second  method  that  may  be  used  to  obtain  solutions  at  the  axis 
is  the  method  of  series  truncations.  The  flow  variables  are  expanded 
into  series  centered  at  the  axis.  In  order  to  obtain  solutions  to  the 
resulting  equations,  it  is  necessary  to  terminate,  or  truncate,  the 
series  after  a  finite  number  of  terms.  Since  the  one-term  truncation 
Is  equivalent  to  the  method  of  local  similarity,  the  series -truncation 
method  may  be  considered  to  be  a  generalization  of  the  method  of  local 
similarity.  Thus  the  two-term  (and  higher)  truncations  can  be  used  to 
test  the  validity  of  the  lccal-similitude  concept.  One  previous  appli¬ 
cation  of  series -truncation  methods  to  the  viscous  blunt -body  problem 
has  been  mane,  by  H-C.  Kao  [21],  and  he  concluded  that  the  local- 
similitude  assumption  was  accurate.  However,  Kao  required  that  the  body 
and  the  shock  be  concentric,  and  thus  the  underlying  assumption  of 
spherical  symmetry  was  left  untested.  The  validity  of  local  similitude 
will  be  re-examined  when  the  series-truncation  analysis  Is  used  to 
obtain  solutions  of  the  simplified  shock-layer  equations  at  the  axis. 

It  will  be  found  that  the  boundaiy  condition  that  is  imposed  on  the 
wall  temperature  plays  a  critical  role  in  this  examination  and,  in 
general,  the  local-similitude  analysis  can  lead  to  substantial  errors. 


In  the  investigation  of  finite-difference  methods,  the  thin-shock- 
layer  analysis  of  H.K.  Cheng  [1,2,3]  serves  as  a  standard  reference. 

In  his  investigation,  the  Havier-Stokes  equations  are  reduced,  under 
the  thin-shock- layer  assumption,  to  a  system  of  parabolic  partial 
differential  equations.  These  equations  are  then  solved  by  an  implicit 
finite -difference  procedure.  The  equations  to  be  investigated  in  the 
subsequent  chapters  contain  several  factors  which  will  make  the  finite- 
difference  analysis  more  difficult  than  the  thin-shock-layer  analysis. 
These  complications  are  associated  with  the  description  of  the  inviscid 
pressure  field  and  of  the  shock  location  when  the  thin-shock-layer 
assumptions  are  removed.  These  factors  have  been  investigated  by  Davis 
and  Chyu  [ll]  for  a  constant -density  shock  layer,  and  they  found  that 
it  was  necessary  to  retain  some  features  of  the  thin-shock-layer  model. 
The  extension  of  the  analysis  of  Davis  and  Chyu  to  flows  of  variable 
density  will  be  considered  here,  as  well  as  two  other  methods  used  on 
implicit  finite -difference  procedures. 

In  the  investigation  that  has  been  outlined  above,  the  problem  is 
formulated  as  a  direct  one  --  that  is,  the  body  shape  is  given  and  the 
shock  shape  is  calculated.  The  alternative  method  of  computing  the 
body  shape  from  a  known  shock  shape  has  been  the  basis  of  numerous 
studies  of  the  blunt-body  problem.  Ihis  inverse  blurt-body  problem 
generally  leads  to  a  considerable  reduction  of  the  computational  diffi¬ 
culties  but  requires  several  computations  if  a  particular  body  shape  Is 
to  be  obtained-  It  is  clearly  shown  in  the  following  chapters  that  the 
shapes  of  the  shock  and  the  body,  relative  to  one  another,  are  quite 
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dependent  upon  the  boundary  conditions  that  are  imposed  at  the  body 
surface  for  the  flow  of  a  viscous,  heat-conducting  fluid.  This  feature 
is  very  important  to  the  use  of  the  series-truncation  or  local-similarity 
methods  --  although  the  inverse  blunt-body  problem  would  reduce  the 
computational  difficulties,  it  would  do  so  at  the  cost  of  leaving  some 
ambiguity  as  to  what  problem  has  been  solved.  Hence,  the  direct  blunt- 
body  problem,  despite  the  inherent  computational  difficulties  that  are 
associated  with  it,  is  considered  to  be  more  appropriate  for  the  present 
investigation  than  the  inverse  problem  is. 

In  the  above  discussion,  only  those  references  of  immediate  concern 
to  the  analysis  of  the  following  chapters  have  been  cited.  For  a  more 
complete  account  of  the  work  done  on  the  blunt-body  problem,  the  reader 
is  referred  to  the  reviews  by  Probstein  [27],  Van  Dyke  [33],  and  Cheng 
[2,3]. 
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Chapter  II 

FORMULATION  OF  THE  PROBLEM 

The  problem  considered  in  this  work  is  that  of  the  flow  about  a 
blunt  body,  which  ia  assumed  to  be  everywhere  convex  and  to  have  a  con¬ 
tinuous  slope.  The  flow  may  be  either  axisymmetric  or  plane  symmetric 
slthough  only  axisymmetric  flows  are  considered  in  the  application  of 
the  equations  to  specific  problems.  A  simplified  form  of  the  Navier- 
Stokes  equations  is  presented  that  is  uniformly  valid  throughout  the 
entire  shock  layer.  The  nsture  of  these  equations  is  examined,  and,  for 
comparison,  some  of  the  features  of  the  equations  used  in  several  thin- 
shock-layer  Investigations  are  discussed. 

A.  Coordinate  System 

A  blunt  body  lies  in  a  flow  field  that  has  a  constant  free-stream 
velocity,  U*  ,  parallel  to  the  body  axis  and  a  density  and  temperature 
given  by  p*  and  T*  respectively.  The  body  surface  is  located  a  dis¬ 
tance  r*  from  the  axis  of  symmetry  and  is  inclined  at  an  angle  ^  -  0 
to  the  axis  ss  shown  in  figure  2.1.  The  body  has  a  longitudinal  curv- 
sture  given  by  K*  and  a  nose  radius  of  a*  .  The  coordinates,  s*  and 
n*  ,  measure  the  arc  length  along  the  surface  of  the  body  snd  the  dis¬ 
tance  normal  to  it.  The  velocity  of  the  fluid  is  resolved  into  compo¬ 
nents,  u*  and  v*  ,  parallel  to  s*  snd  n*  respectively.  A  shock 
of  zero  thickness  is  located  at  a  distance  A*  from  the  body.  The  angle 
between  the  shock  and  a  plane  normal  to  the  axis  is  denoted  by  0. 
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B,  Non-dimensional  Variables 

Dimensionless  variables  are  introduced  by  referring  all  lengths  to 

2 

a*  ,  velocities  to  U*  »  «.ne  density  to  p*  »  pressure  to  o*U*  , 

CD  CD  CD 

2 

temperature  to  U*  /C*,  and  viscosity  to  --its  value  at  the  refer- 

2  1 

ence  temperature  U*  /C*  .  With  such  a  choice  of  reference  values,  the 

CD  p 

variables  all  remain  bounded  in  the  hypersonic  limit,  M  ->»  .  In 

00 

addition,  the  normal  coordinate  and  normal  velocity  are  stretched  by  a 
factor  1/e  where 


e  ■ 


(2.1) 


i 

Starred  symbols  denote  dimensional  quantities  and  unstarred  sym¬ 
bols  denote  dimensionless  quantities. 
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The  normal  coordinate*  n  ,  and  the  normal  velocity,  v  ,  are  iden¬ 
tical  in  form  to  those  used  in  boundary- layer  theory.  This  form  was 
chosen  simply  as  a  convenience  for  the  numerical  solution  of  the  equa¬ 
tions.  Since  the  viscous  effects  are  confined  to  a  region  of  order  e 
near  the  body,  the  use  of  the  stretched, boundary-layer  variables  guaran¬ 
tees  that  the  region  of  large  velocity  gradients  does  not  become  exces¬ 
sively  thin  as  e  becomes  small.  However,  tne  use  of  this  form  of  the 
variables  is  not  essential  to  the  method  of  solution  developed  here. 


C.  Shock-layer  Equations 

The  Navier-Stokes  equations  can  be  written  in  the  orthogonal, 
curvilinear  coordinate  system  described  above  (see,  for  example,  ref¬ 
erence  32  ).  The  quantity  e  ,  defined  by  equation  (2.1;  and  referred 
to  as  the  viscous,  hypersonic  similarity  parameter,  appears  in  the  momen¬ 
tum  and  energy  equations  when  the  dimensionless  variables  of  section  B 
are  Introduced.  In  the  hypersonic  limit,  it  is  the  only  similarity  par¬ 
ameter  that  appears  in  the  formulation  of  the  problem.  The  equation  for 
the  similarity  parameter  may  be  rewritten  in  terms  of  the  free-stream 
Mach  number  and  Reynolds  number  by  using  the  relation 


U* 

<0 

C* 

P 


(Y-l)  M  T* 

'  00  CD 


If  a  power  law  is  used  for  the  viscosity,  y*  a  T*^  ,  then 


e  -  (y-l) 


W2  « 


In  this  form,  it  can  be  seen  that  e  is  a  measure  of  the  ratio  of  the 
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mean  free  path  in  the  shock  layer  to  the  thickness  of  the  shock  layer 
[18,  p.378].  Thus,  e  aerves  as  a  measure  of  the  degree  of  rarefaction 
of  the  flow  in  the  shock  layer.  Further,  the  quantity  e  is  the  pertur¬ 
bation  parameter  used  to  systematically  expand  the  Navier-Stokes  equa¬ 
tions  into  the  equations  of  boundary- layer  theory  [  32  J.  Thus  the 
second-order  boundary -layer  theory  represents  an  extension  of  the  con¬ 
tinuum  flow  equations  toward  the  regime  of  rarified  flow.  Adequate 
numerical  techniques  for  the  Ljlut  '.on  of  th<;  first-  and  second-order 
boundary- layer  equations  exist  (e.g. ,  Flugge-Lotz  and  Blottner  [  14], 
Davis  and  Fliigge-Lotz  [12],  Clutter  and  Smith  [  5  ]),  but  there  is 
still  some  difficulty  associated  with  the  application  of  boundary-layer 
theory.  In  particular,  it  requires  the  solution  of  a  series  of  inter¬ 
related  problems,  two  of  which  (the  basic  inviscid  flow  and  the  invis- 
ctd  flow  due  to  the  displacement  thickness)  involve  solving  elliptic 
partial  differential  equations  in  regions  in  which  the  locations  of  the 
boundaries  are  unknown. 

Numerous  methods  have  been  developed  to  solve  this  inviscid  flow 
problem.  Perry  and  Pasiuk  requested  solutions  to  a  few  standard  prob¬ 
lems  from  various  agencies  which  have  developed  numerical  procedures 
based  on  some  of  these  methods.  Thfs  survey  included  solutions  obtained 
from  both  direct  and  inverse  methods  and  from  integral-relation  methods, 
finite-difference  methods,  series-truncation  procedures,  and  a  time- 
dependent  method.  A  brief  comparison  of  the  results  of  this  survey  is 
presented  in  reference  25  .  Several  of  ''e  solutions  predicted  surface 
pressures  that  are  in  excellent  agreem.-'  with  experimental  data;  but, 
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in  general,  considerable  disparity  exists  in  the  results. 

Thus  it  appears  that  accurate  results  to  the  inviscid  flow  problem 
can  be  obtained  with  sufficient  numerical  skill,  but  these  results  must 
then  be  used  in  an  entirely  different  numerical  procedure  to  obtain  the 
first-order  boundary-layer  results.  Then,  before  all  the  second-order 
boundary- layer  effects  can  be  computed,  the  inviscid  flow  due  to  the  dis¬ 
placement  thickness  effect  of  the  boundary  layer  must  be  obtained  either 
by  solving  the  inviscid  flow  problem  again  or  by  using  some  approximation 
(c.f.  ref.  12  )•  Thus  the  application  of  higher-order  boundary-layer 
theory  is  complicated  by  the  necessity  of  solving  several  different  sets 
of  equations  which  require  different  methods  of  solution. 

As  an  alternative,  Davis  and  Fliigge-Lotz  [  12  ]  proposed  representing 
the  entire  shock  layer  by  one  set  of  equations,  a  simplified  form  of  the 
Navier-Stokes  equations  that  would  be  uniformly  valid  throughout  the 
shock  layer  and  would  contain  all  the  viscous  effects  that  are  contained 
in  the  boundary- layer  equations.  Although  a  solution  to  this  set  of 
equations  may  require  the  use  of  elaborate  numerical  techniques,  only 
one  application  of  the  numerical  methods  is  needed.  This  paper  repre¬ 
sents  an  investigation  of  this  alternative  procedure  proposed  by  Davis 
and  FlUgge-Lotz.  To  obtain  the  desired  equations,  only  those  terms  of 
the  Navier-Stokes  equations  that  contribute  to  che  first-  or  second- 
order  boundary- layer  equations  or  contribute  to  the  first-  or  second- 
order  outer,  inviscid  flow  are  retained.  Thus  the  modified  equations 
are  uniformly  valid  to  order  e  throughout  the  shock  layer.  The 
shock- layer  equations  that  are  obtained  in  this  manner  are 


Continuity 


[(r  +  en  sin  0)^  pul  +  [(1  +  <  €*0  (r  +  en  sin  0)^  pv)  =  0 
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The  subscripts  s  and  n  denote  partial  differentiation  with 
respect  to  the  s  and  n  coordinates,  respectively.  The  quantity  j 
is  equal  to  0  for  plane  flow  and  equal  to  1  for  axlsymmetrlc  flow. 


In  equation  (2.5)  the  Prandtl  number,  (7  ,  has  been  assumed  to  be  a 
constant,  but  this  assumption  is  not  essential  to  the  use  of  the  methods 
of  analysis  investigated  in  later  chapters. 

The  thermodynamic  state  relation  is  taken  to  be  as  simple  as  possible 
in  order  to  facilitate  computations.  Thus  a  perfect  gas  with  constant 
specific  heats  is  assumed.  Hence, 

State 


This  assumption  represents  a  limitation  on  the  applicability  of  the  results 

of  this  work.  For  some  flow  conditions,  it  is  necessary  to  consider 
the  real  gas  effects  arising  from  the  molecular  structure  of  the  gas  or 
from  chemical  reactions  among  the  various  components  of  the  gas.  In 
addition,  the  reaction  times  for  such  molecular  processes  become  signifi¬ 
cant  for  some  problems,  and  in  such  a  case  it  is  necessary  to  account  for 
the  fact  that  the  chemical  state  of  the  gas  may  not  be  in  equilibrium. 
However,  it  is  desired  in  this  work  to  concentrate  on  the  difficulties 
arising  solely  from  the  hydrodynamic  aspects  of  the  flow, and  therefore 
all  complications  due  to  the  chemical  nature  of  the  gas  have  been  ignored. 
To  the  above  equations  must  be  added  an  equation  which  gives  the  vis¬ 
cosity  as  a  function  of  the  temperature: 

p,  *  p,(T)  (2.  7) 
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Although  equations  (2.2)-(2.7)  have  been  chosen  simply  as  an  attempt 
to  have  a  somewhat  simpler  formulation  of  the  problem  than  is  afforded 
by  the  boundary-layer  theory,  there  Is  one  aspect  of  the  flow  In  which 
the  current  formulation  of  the  problem  may  be  superior  to  that  of 
boundary- layer  theory.  Recently  Weinbaum  and  Garvine  140]  have  Investi¬ 
gated  the  nature  of  the  pressure  field  In  compressible  boundary -layer 
theory.  They  assert  that  under  certain  flow  conditions  there  are  three 
distinct  regions  requiring  description  In  a  viscous  flow  problem  with 
large  Reynolds  number.  Near  a  solid  boundary,  the  usual  boundary  layer 
exists  where  the  normal  pressure  gradient  Is  negligible  and  the  viscous 

effects  predominate.  Far  from  the  boundary,  there  Is  the  usual  invlscid 

1 

region.  However,  they  state  that  when  Mg  k  >  1  there  Is  an  Inter¬ 
mediate  region  where  the  effects  of  viscosity  are  considerable  and  where 
the  normal  pressure  gradient  cannot  be  neglected.  Thus,  they  propose 
the  existence  of  a  transition  region  between  the  boundary  layer  and  the 
Invlscid  flow.  Under  such  conditions,  the  matching  of  the  boundary- 
layer  solution  and  the  invlscid  solution  as  done  In  the  development  of 

3 

higher-order  boundary -layer  theory  may  not  be  valid. 

It  should  be  noted  that  this  proposed  Intermediate  region  between 

the  boundary  layer  and  the  invlscid  flow  would  not  occur  in  the  subsonic 

portion  of  the  flow  around  a  blunt  body.  However,  Davis  and  Flugge-Lotz 

1 

The  subscript  e  refers  to  a  streamline  at  the  outer  edge  of  the 

viscous  layer. 

2 

It  has  also  been  suggested  [38]  that  this  Is  a  "strong-interaction" 
phenomenon,  which  properly  lies  outside  the  scope  of  boundary-layer 
theory. 


I 


[  12]  report  that  in  the  flow  around  a  sphere  the  normal  velocity  found 
from  the  second-order  boundary- layer  theory  grows  extremely  rapidly  (com¬ 
pared  with  the  first-order  normal  velocity)  after  the  invlscid  sonic  line 
is  passed.  Pavis  attributed  this  to  a  breakdown  of  boundary- layer  theory 
as  the  boundary-layer  separation  point  is  approached.  Fannel'dp  and 

FlUgge-Lotz  [  13  ]  report  a  similar  result  for  flow  around  a  circular 
l 

cylinder.  However,  Fannel'dp  suggests  that  the  result  is  not  a  manifes¬ 
tation  of  the  physical  singularity  that  occurs  at  separation  but  that 
it  is  a  singularity  due  to  a  limitation  of  the  mathematical  theory. 

Further,  he  suggests  that  the  limitation  is  related  to  an  Inadequate 
description  of  the  streamline  curvature. 

2 

The  normal  pressure  gradient  is  proportional  to  p  k  u  in 

second-order  boundary -layer  theory.  This  represents  the  centrifugal 

force  which  occurs  in  a  flow  field  whose  streamlines  have  the  same  local 

curvature  as  the  body.  Fannel'dp  points  to  this  similarity  to  the  thin- 

shock  layer  model  and  suggests  hat  the  difficulty  could  be  analogous  to 

the  occurrence  of  the  zero-pressure  point  [18,  p.82]  in  thin-shock- layer 

theory.  This  is  not  inconsistent  with  the  analysis  of  Weinbaum  and  Garvine 

since  the  interaction  of  the  normal  pressure  gradient  and  local  variations 

in  the  streamline  curvature  (represented  by  the  terms  proportional  to 

u  v  and  v  v  in  the  normal  momentum  equation)  plays  an  essential  role 
s  n 

in  their  analysis.  The  set  of  equations  (2. 2) - (2. 7)  possesses  a  complete 
description  of  the  Interaction  of  the  pressure  and  the  local  streamline 
curvature  and  therefore  should  not  encounter  this  difficulty. 

1  For  non-circular  bodies  for  which  K<  1  downstream  of  the  sonic 
line,  no  difficulties  were  reported  in  either  reference  12  or  13. 
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Weinbaum  and  Garvine  also  adopt  the  equations  proposed  by  Davis  and 
Flugge-Lotz,  and  they  propose  a  method  of  solution.  Weinbaum  [39  ]  has 
appltc  this  method  to  solve  for  the  flow  in  a  laminar  wake.  The  use  of 
.is  method  to  solve  the  blunt-body  problem  is  investigated  in  Chapter 
IV. 


D.  Characteristics  of  the  Equations 

It  is  interesting  to  examine  the  characteristics  of  the  set  of  equa¬ 
tions  (2. 2) -(2. 7).  The  characteristics  fall  into  two  categories.  Most 

i 

of  the  characteristics  coincide  and  are  described  by 

a  *  constant  (2.8) 

Although  the  classifications  of  hyperbolic,  parabolic,  and  elliptic  are 
not  strictly  applicable  to  such  higher-order  equations,  these  character¬ 
istics  may  be  thought  of  as  being  parabolic  in  nature.  This  set  of  char¬ 
acteristics  can  be  traced  to  the  tangential  momentum  equation  and  the 
energy  equation.  Thus  these  two  equations,  which  describe  the  behavior 
of  u  and  T  (u  and  T  occur  only  in  these  two  equations),  are 
still  essentially  parabolic  just  as  they  are  in  boundary -layer  theory; 
but  now  the  pressure  gradient  is  not  a  known  function  but  is  determined 
together  with  v  from  the  interaction  of  the  continuity  and  normal 
momentum  equations.  The  remaining  characteristics  are  described  by 

£  dn  _  ev  (  /Fl  _  ev  ,  £  (2.9) 

1  +  e/cnds  u  — p  u  u 

- 1 - 

The  equation  describing  the  characteristic  surfaces  for  a  set  of 
nonlinear  partial  differential  equations  may  be  obtained  from  reference 
26  ,  page  32. 
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That  is,  there  are  two  additional  characteristics  whose  slopes  differ 
from  that  of  a  streamline  by  an  amount  •  Thus  the  behavior  of 

the  normal  momentum  and  continuity  equations  is  essentially  hyperbolic. 
Despite  the  difficulty  of  dealing  with  a  system  of  such  mixed  character, 
it  would  appear  that  the  equations  should  be  well  suited  for  an  initial- 
value  problem. 

Since  we  know  that  the  subsonic,  inviscid  flow  region  of  the  shock 
layer  is  elliptic  in  character,  it  is  necessary  to  inquire  as  to  how  or 
whether  this  "parabolic-hyperbolic"  set  of  equations  can  adequately 
describe  the  flow  in  a  region  in  which  some  upstream  influence  should 
occur.  Weinbaum  [40  ]  has  reported  that  these  equations  are  stable  as 
an  initial-value  problem  when  applied  to  laminar-wake  flows.  However, 
even  though  the  modified  equations  are  not  of  an  elliptic  nature  and  thus 
may  be  numerically  stable  in  an  initial-value  problem,  one  would  not 
expect  the  inherent  elliptic  character  of  the  blunt -body  problem  to  be 
absent.  The  manner  in  which  an  upstream  influence  can  occur  is  not  read¬ 
ily  apparent, and  the  answer  to  this  must  be  deferred  until  the  methods  of 
solution  are  investigated. 
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H.  K.  Cheng  has  extensively  Investigated  this  flow  model  In  a  series 

of  papers  [  1,  2,  3  ].  The  thln-shock-layer  equations  are  an  asymp- 

p® 

totlc  description  of  the  shock  layer  for  small  density  ratio,  -  . 

pah 

Cheng  uses  the  first-order  equations  that  result  from  an  order  of  mag- 

nltude  analysis  for  small  -  .  These  equations  are  parabolic  In  nature 

psh 

(all  of  the  characteristics  coincide)  and,  In  fact,  are  Identical  to  the 

boundary- layer  equations  except  that  the  normal  pressure  gradient  Is 

2 

proportional  to  p  <u  .  This  enables  Cheng  to  solve  the  equations  with 
an  Implicit  finite-difference  scheme  which  starts  at  the  axis  of  symmetry 
and  marches  downstream.  This  property  of  the  equations  Is  quite  desirable 
since  Implicit  finite-difference  methods  have  been  well  developed  and 
lead  to  an  accurate  solution  over  as  large  a  region  as  is  desired.  Cheng 
docs  not  assume  an  Infinitesimally  thin  shock  but  uses  a  modified  set 
of  the  Rank lne-Hugon lot  conditions  to  account  for  transport  effects  at 
the  shock.  These  relations  are  similar  to  those  proposed  by  Sedov, 
Michallova  and  Chernyl  [  29  ].  One  important  simplification  In  the  bound¬ 
ary  conditions  results  from  the  use  of  the  thln-shock-layer  concept.  To 
the  order  of  the  approximation  considered  by  Cheng,  there  is  no  distinc¬ 
tion  between  the  location  of  the  shock  Interface  and  the  body  surface. 

This  appears  in  most  of  the  analyses  based  on  the  thin  shock  la>?r  as  an 
assumption  that  the  body  and  shock  are  concentric  (0  *  0)  .  This  sim¬ 
plifies  the  analysis  considerably  but,  of  course,  influences  the  accuracy 
of  the  results.  Despite  the  simple  nature  of  this  flow  model,  Cheng 
obtains  reasonable  agreement  with  the  results  of  other  Investigations. 
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1 

However,  since  several  terms  of  0(e)  ,  including  the  slip  and  temper¬ 

ature  Jump  effects,  are  missing  from  his  formulation,  there  is  some  doubt 
as  to  the  accuracy  of  the  solutions  at  the  low  Reynolds  numbers  that 
were  considered.  Cheng  avoids  any  potential  difficulty  due  to  the 
occurrence  of  a  zero-pressure  point  by  considering  only  flow  around 
such  bodies  as  paraboloids  and  hyperboloids.  However,  as  discussed  in 
the  previous  section,  it  would  be  desirable  to  have  a  more  complete 
description  of  the  pressure  field. 

Other  investigations  have  retained  more  terms  of  the  equations  than 
Cheng  did.  This  provides  a  more  accurate  description  of  the  shock  layer 
at  low  Reynolds  numbers  but  increases  the  difficulty  of  obtaining  accurate 
solutions  since  the  equations  are  no  longer  parabolic.  The  method  of 
solution  usually  employed  in  such  Investigations  is  to  reduce  the  equa¬ 
tions  to  ordinary  differential  equations  by  assuming  the  functional 
dependence  of  the  variables  with  the  s-coordinate.  Examples  of  such 
analyses  are  those  of  Ho  and  Frobsteln  [19  J,  Shih  and  Krupp  [30  J,  and 
L.  Goldberg  [17  j. 

The  investigation  of  Ho  and  Probstein  removed  the  constant  density 
assumption  which  had  characterized  many  prior  investigations,  and 
many  subsequent  investigations  have  followed  the  same  basic  approach. 

Ho  and  Probstein  emphasize  that  the  thln-shock-layer  approximation  is  not 
essential  to  their  investigation  but  is  used  only  to  make  the  numerical 
analysis  simpler.  .  This  is  true  of  most  of  the  thin-shock-layer  investi- 

- 1 - 

A  function  f(|)  is  0(|)--of  order  $--if  the  limit  of  f(|)/| 
exists  as  |  -»  0  . 
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gations  since  the  equations  are  generally  not  reduced  to  the  parabolic 
system  used  by  Cheng.  Even  though  the  use  of  the  thin-shock-layer  equa¬ 
tions  is  not  essential  to  these  investigations,  the  consistent  assumption 
that  the  body  and  shock  are  concentric  is  very  important.  These  investi¬ 
gations  assume  that  the  solution  has  a  locally  similar  character  in  the 
stagnation  region.  This  assumption  is  based  on  the  existence  of  spheri¬ 
cal  or  cylindrical  symmetry  when  the  body  and  shock  are  concentric.  H.  C. 
Kao  [21  ]  has  tested  the  validity  of  the  local  similarity  concept  by 
solving  the  Navier-Stokes  equations  by  the  method  of  series  truncations. 
Kao  concluded  that  the  similarity  solutions  were  quite  accurate.  How¬ 
ever,  it  should  be  noted  that  Kao  also  assumed  that  the  body  and  shock 
were  concentric.  Thus  the  underlying  assumption  of  spherical  symmetry 
was  left  untested.  This  question  is  examined  in  Chapter  111  when  the 
method  of  series  truncation  is  considered. 

The  investigetion  of  Goldberg  [  17  ]  is  similar  to  that  of  Ho  and 
Probstein  but  uses  a  somewhat  more  elaborate  flow  model.  In  particular, 
Goldberg  modifies  the  shock  conditions  in  a  manner  similar  to  that  pro¬ 
posed  by  Sedov,  et  al.  However,  the  basic  assumption  that  the  shock 
interface  and  the  body  are  concentric  remains. 

W.  C.  L.  Shih  and  R.  S.  Krupp  [  30  ]  have  also  used  the  thin-shock- 
laycr  assumption  to  simplify  the  Navier-Stokes  equations  but  have  not 
assumed  that  the  shock  slope  is  known.  In  addition  the  method  of  analy¬ 
sis  is  somewhat  modified  in  order  to  obtain  the  values  of  the  variables 
at  a  specified  number  of  positions  downstream  of  the  axis  of  symmetry. 
Their  analysis  also  provides  for  the  systematic  improvement  of  the  degree 
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of  approximation.  In  each  of  the  above  investigations,  except  that  of 
H.  K.  Cheng,  the  terms  proportional  to  puv  and  pvv  have  been 

8  O 

retained  in  the  normal  momentum  equation.  Thus  a  complete  description 
of  the  inviscid  pressure  field  is  afforded.  However,  these  investiga¬ 
tions  retain  a  different  set  of  viscous  terms  than  has  been  retained  in 

equations  (2. 2)- (2. 7).  In  particular  the  term  (u  v  )  is  retained  in 

n  n 

the  normal  momentum  equation.  This  term  has  considerable  Influence  on 
the  character  of  the  equations.  Shih  and  Krupp  state  that  such  a  system 
of  equations  is  parabolic.  Davis  and  Flugge-Lotz  had  also  suggested  the 
possibility  of  reducing  the  equations  to  parabolic  form  by  Including  the 
vn)n  tern  in  the  normal  momentum  equation.  However,  the  inclusion 
of  this  term  results  in  the  streamlines  being  characteristics, and  it  is 
not  clear  that  the  equations  should  be  classified  as  parabolic  in  such  a 
case.  Since  the  streamlines  become  parallel  to  the  body  as  n  -* 0,  the 
body  surface  will  be  a  characteristic  surface  for  such  a  set  of  equations. 
Thus  by  the  definition  of  characteristics,  one  cannot  solve  for  the  high¬ 
est-ordered  derivatives  normal  to  the  body,  given  all  other  flow  quanti¬ 
ties  on  the  body,  unless  the  specified  boundary  conditions  satisfy  a  com¬ 
patibility  condition  along  the  boundary.  If  the  boundary  conditions  do 

have  this  special  form,  then  there  is  no  guarantee  of  the  uniqueness  of 
l 

the  solution.  This  fact  is  very  Important  to  any  method  of  analysis  that 
reduces  the  problem  to  one  of  solving  ordinary  differential  equations 
along  lines  normal  to  the  body.  (This  point  will  be  illustrated  in 

Chapter  III.)  Shih  and  Krupp  modify  their  equations  by  neglecting 
- 1 - 

Discussions  of  these  points  may  be  found  in  reference  16  , 
especially  sections  4.1,  2.1,  3.1  and  3,4. 
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viscous  forces  normal  to  streamlines  in  order  to  avoid  au  instability 

which  occurs  near  the  body  surface.  This  modification  does  not  remove 

the  (u  v  )  term,  but  does  alter  its  special  role  described  above, 
n  n 

Kao  also  reports  an  instability  near  the  body  surface, and  he  consequently 
uses  a  modified  set  of  equations  near  the  surface.  Ho  and  Probsteln  use 
a  knowledge  of  the  general  properties  of  the  variables  near  the  wall  to 
eliminate  any  numerical  difficulty  with  the  singularity.  It  seems  likely 
that  such  difficulties  are  related  to  the  attempted  use  of  a  character¬ 
istic  surface  as  a  boundary  surface. 

In  addition,  Shih  and  Krupp  have  computed  an  alternate  solution  to 
an  example  solved  by  Ho  and  Probsteln,  verifying  the  non-uniqueness  which 
may  occur  from  such  a  formulation  of  the  problem. 

F.  Boundary  Conditions  at  the  Body  Surface 

Hayea  and  Probsteln  [  18  ]  state  that  the  effects  of  slip  and  tem¬ 
perature  jump  on  heat  transfer  and  wall  shear  are  negligible  provided 
that  the  wall  temperature  is  low  ,  even  if  the  shock  layer  cannot  be  con¬ 
sidered  a  continuum.  All  the  investigations  cited  in  the  preceding 
section  have  used  no-slip  boundary  conditions  for  this  reason.  In  this 
present  work,  however,  the  wall  boundary  conditions  will  take  the  slip 
and  temperature -jump  effects  into  consideration  since  they  are  0(e)  in 
the  boundary  layer  and  therefore  should  be  included  to  be  consistent  with 
the  retention  of  terms  of  0(e)  in  the  partial  differential  equation!). 
Thus  the  examples  to  be  calculated  will  not  be  restricted  to  very  low 
wall  temperatures.  Davis  and  Fltigge-Lotz,  in  their  study  of  the  second- 
order  boundary -layer,  have  found  that  the  effect  on  the  stagnation-point 
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heat  transfer  is  not  negligible  for  moderately  cold  walls. 
The  wall  boundary  conditions  are 


u(s,0) 


e  mT_(s,0)] 

P(«,0) 


0 


(2.10) 


T(s,0) 


Tb(8)  +  e 


IAtTC«,0>  J 

p(s,0) 


- - - - 1 

T(s ,0) 


dn 


0 


(2.11) 


v(s,0)  »  0 


(2.12) 


a^  and  c^  are  dimensionless  constants  which  are  0(1).  These  condi¬ 
tions  are  seen  to  be  those  given  by  Van  Dyke  [32  ],  except  that  a  term 

jyr 

proportional  to  §=■  has  been  omitted  from  equation  (2.10)  since  It  Is 
08 

2 

0(e  )  In  the  boundary  layer.  For  comparison  with  the  results  of  the 
second-order  boundary- layer  Investigation  of  Davis  and  Flugge-Lotz,  the 
values  of  a^  and  c^  are  taken  to  be  and  g^( ^ ^  ^  respec¬ 

tively.  Equation  (2,11)  may  be  replaced  with  a  suitable  condition  on 
the  heat  transfer.  In  Chapter  III,  several  examples  will  be  considered 
In  which  an  adiabatic  wall  Is  specified.  In  this  case  equation  (2.11) 

Is  replaced  with 


T  (s,0)  +  a  u(s,0)  u  (s,0)  •  0  (2.11a) 

n  n 

G.  Shock  Conditions 

As  stated  earlier,  the  shock  Is  assumed  to  be  vanishingly  thin, and 
therefore  the  Ranklne-Hugonlot  conditions  may  be  applied.  These  condi¬ 
tions  may  be  written  In  a  form  In  which  the  values  of  the  variables 
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Immediately  behind  the  shock  are  given  explicitly: 


u(8,A)  -  cos(0-0)  sin  0  +  ■  *.2  sin(9-0)  cos  0 

P  ( 8 » A) 


(2.13) 


ev(s,A)  ■  sin  0  sin(9-0) - —  cos  0  cos(0-0) 

p(s,A) 


P(s, 


p(s,A) 


-  2  r  2  Y-i  n 

*  *  cos 

^  00  _ 


—  -  + 
A\  Y+l 


(Y+l)  M2  cos2  0 

00 


T(s,A) 


Y-l  ,  2 

2 

2 

y  cos  a  1 

Y+l  (Y+l)  M2  cos2  0 

00 

Y+l 

Y-l  '  2  M2 

00 

(2.14) 

(2.15) 

(2.16) 

(2.17) 


Equations  (2. 15)-(2. 17)  are  not  independent,  of  course,  but  are  related 
through  the  state  equation  (2.6).  The  factors  of  cos(0-0)  and 
sin(0-0)  in  equations  (2.13)  and  (2.14)  represent  a  rotation  of  axes  to 
change  from  velocity  components  tangential  and  normal  to  the  shock  to 
components  tangential  and  normal  to  the  body  surface. 

The  assumption  of  a  thin  shock  is  consistent  with  the  rest  of  the 

2 

mathematical  model  proposed  here  as  the  shock  thickness  is  0(e  ) 

(18  ,  32]  and  thus  need  not  be  considered  here. 

It  is  advantageous  at  times  to  replace  one  of  the  Rankine-Hugoniot 
conditions  with  an  equation  that  expresses  the  overall  balance  of  mass 
flow: 


(r  +  eA  sin  9)^  +  1  *  2^  e  J  p  u  (r  +  en  sin  0)^  dn  (2.18) 
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H.  Additional  Relations 


There  are  several  geometric  relationships  which  are  useful  in  the 
solution  of  the  problem.  The  most  essential  equation  expresses  the 
relation  between  the  shock  position  and  its  slope; 


dA  - 

e  ^  *  (1  +  e  K  A)  tan  (0  -  0) 


(2.19) 


In  addition,  the  body  geometry  gives 


K 


d0 

ds 


(2.20) 


and 

=  cos  0  (2.21) 

The  shear  stress  and  heat  transfer  at  the  body  surface  are  useful 
for  evaluating  solutions  of  the  equations.  The  non-dimensional  shear 
stress  and  heat  transfer  are  given  by 


and 


(2.22) 


(2.23) 


I.  Sumnary 

Equations  (2.2)- (2. 7) ,  (2. 10) - (2. 17)  comprise  the  system  for  which 
we  seek  a  solution.  These  equations  are  uniformly  valid  to  0(e) 


throughout  the  entire  shock  layer  and  therefore  may  be  used  to  describe 
the  flow  when  the  fluid  Is  slightly  rarlfled.  The  equations  should  be 
well  suited  for  an  Initial- value  problem  and  thus  should  avoid  the 
numerical  Instabilities  attendant  to  solving  the  Euler  equations  In  the 
subsonic  portion  of  the  flow.  Having  one  set  of  equations  for  the 
entire  shock  layer  eliminates  the  need  of  solving  a  series  of  Inter¬ 
related  problems  as  Is  done  In  boundary- layer  theory;  but,  of  course, 
the  ability  to  determine  the  Individual  second-order  effects  of  longi¬ 
tudinal  and  transverse  curvature,  vortlclty  Interaction,  etc.  Is  lost. 


Chapter  III 

METHOD  OF  SERIES  TRUNCATION 

In  the  previous  chapter  It  was  noted  that  the  set  of  equations 
describing  the  flow  In  the  shock  layer  were  of  such  a  nature  that  they 
should  be  solvable  as  an  Initial-value  problem.  A  necessary  step  In 
such  a  solution  Is  to  obtain  an  accurate  set  of  Initial  values.  The 
Initial  station  for  this  problem  Is  the  axis  of  symmetry  of  the  flow. 
Since  the  values  of  the  variables  on  the  axis  are  not  known  a  priori, 

It  Is  necessary  to  obtain  a  solution  to  equations  (2. 2) -(2. 7)  that  Is 
valid  In  the  neighborhood  of  the  axis.  In  this  chapter  a  method  Is 
described  by  which  such  a  solution  can  be  found.  The  results  of  the 
computation  of  several  examples  are  presented  In  order  to  evaluate  the 
accuracy  of  the  method  of  solution  and  to  evaluate  the  range  of  validity 
of  the  basic  flow  model.  The  concept  of  local  similarity  as  a  method  of 
analysis  of  blunt -body  flows  Is  reviewed. 

A.  Description  of  Series  Truncations 

The  most  apparent  method  of  obtaining  a  solution  at  the  axis  Is  to 
expand  each  of  the  variables  about  the  axis  Into  a  series  In  the 
a  -  coordinate.  The  coefficients  of  these  series  are  functions  of  the 
normal  coordinate,  n.  When  these  series  are  substituted  Into  equations 
(2. 2) -(2. 7),  the  problem  Is  reduced  to  solving  a  set  of  ordinary  differ¬ 
ential  equations  for  the  coefficients  of  the  series.  Determining  the 
first  coefficient  of  each  of  the  series  gives  the  desired  solution  et 


the  axis.  Determination  of  higher-order  coefficients  extends  the  solu¬ 
tion  to  points  away  from  the  axis.  In  theory,  this  may  be  continued  to 
give  results  as  accurate  as  are  desired  within  the  interval  of  conver¬ 
gence  of  the  basic  series. 

Such  a  method  of  analysis  has  been  widely  employed  for  solving  non¬ 
linear  partial  differential  equations  in  two  independent  variables.  The 
Blasius-series  solution  of  the  boundary- layer  equations  is  a  well-known 
example  of  such  a  procedure.  In  the  Blasius  series,  the  coefficient  of 
the  first  term  of  the  series  may  be  calculated  Independently  of  the 
other  coefficients.  The  coefficient  of  any  term  of  the  series  beyond 
the  first  depends  only  on  the  lower-ordered  terms.  Thus  the  coefficients 
may  be  calculated  successively,  starting  with  the  first,  until  the 
desired  degree  of  accuracy  is  attained  throughout  the  region  of  Interest. 

In  the  present  problem,  however,  the  ordinsry  differential  equations 
for  the  coefficients  of  the  series  may  not  be  solved  so  simply.  Due  to 
the  nature  of  the  problem,  the  first  coefficients  of  the  series  are  not 
independent  of  the  higher -order  coefficients.  As  a  consequence,  any 
finite  number  of  the  ordinary  differential  equations  contain  more  unknown 
variables  than  there  are  equations.  Thus,  in  order  to  solve  the  equa¬ 
tions,  it  is  necessary  to  Introduce  some  additional  approximations  for 
the  higher  coefficients.  The  procedure  used  to  handle  this  situation  is 
known  as  the  method  of  series  truncation  and  is  described  below. 

The  variables  are  expanded  into  the  following  series: 

3 

u(s,n)  “  u^(n)  sin  8  +  u^(n)  sin  8  +  ...  (3.1a) 

v(s,n)  «  -  cos  8  [v  (n)  +  v2(n)  sin2  8  +  . . . J  (3.1b) 
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p(s,n)  ■  pQ(n)  +  P2(n)  8in2  $  +  F^(n)  sin^  0  +  ...  (3.1c) 

T(8,n)  «■  T  (n)  +  T_(n)  sin^  0  +  ...  (3.  Id) 

o  i. 

2 

p(s,n)  -  pQ(n)  +  p2(n)  sin  0  +  ...  (3.1e) 

where  the  relation  between  0  and  8  depends  on  the  body  shape.  The 
particular  form  of  the  expansions  is,  in  general,  quite  important  for 
the  ultimate  success  of  the  method  of  series  truncation.  Expansions 
quite  similar  to  those  of  equations  (3.1)  have  been  used  in  previous 
investigations,  and  they  are  reported  to  have  yielded  favorable  results 
for  flow  around  spherical  bodies  [21  ].  For  other  body  shapes,  if  one 
expects  accurate  results  over  a  significant  distance  in  the  s-dlrectlon, 
the  expansion  should  probably  be  chosen  to  take  advantage  of  the  indi¬ 
vidual  characteristics  of  that  particular  problem.  For  example,  various 
types  of  flows  around  a  paraboloid  of  revolution  have  been  analysed  with 
excellent  results  by  transforming  to  parabolic  coordinates  and  expanding 
into  a  suitable  series  in  these  new  coordinates  [10,31,36,4  ].  Although 
no  answer  can  be  given  to  the  question  of  what  form  of  expansions  pro¬ 
vides  the  best  results  for  an  arbitrary  body  shape,  the  expansions  (3.1) 
should  yield  reasonable  results  at  the  axis  of  symmetry  for  any  body. 

The  use  of  0  rather  than  a  in  equations  (3.1)  is  recommended 
since  the  expansion  of  the  governing  equations  is  expected  to  be  some¬ 
what  simpler  with  the  use  of  0.  The  arc  length  s  does  not  appear 
explicitly  in  the  equations  (except  in  which  can  be  written  as 

K  ).  On  the  other  hand,  0  appears  in  both  the  partial  differen¬ 

tial  equations  and  the  shock  boundary  conditions.  Thus  an  expansion 
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in  powers  of  sin(s)  would  require  that  8  be  expanded  in  powers  of 

sin(s).  The  resulting  equations  would  then  be  correspondingly  more 

complex  (the  shock  conditions  become  particularly  difficult  to  handle). 

In  addition,  the  geometric  description  of  the  body  is  frequently  simpler 

in  terms  of  6  ,  which  again  simplifies  the  expansion  of  the  equations. 

(As  an  example,  consider  the  parabola.  In  terms  of  0  ,  we  have 

3 

r  *  tan  8  and  k  *  cos  0.  In  terms  of  the  arc  length  s  ,  we 
have  only  implicit  relations  such  as  s  -  j|r  Vr^  +  1  +  log(r  ’TI,]., 
This  distinction  disappears  for  a  spherical  body  since  8  ■  s  for  the 
sphere. 

The  sphere  has  been  extensively  treated  in  the  literature  despite 
certain  unpleasant  aspects  of  the  flow  around  such  bodies  (f]ow  separa¬ 
tion,  trailing  wake,  etc.).  Since  this  choice  of  body  shape  provides 
numerous  opportunities  for  comparison  with  the  results  of  other  investi¬ 
gations,  the  rest  of  the  work  in  this  chapter  will  be  directed  toward 

l 

the  analysis  of  the  flow  around  spherical  bodies.  This  choice  is  not 
essential  to  the  method  of  analysis,  however,  since  the  expansions  which 
follow  could  have  been  made  for  some  other  choice  of  body  shape. 

For  a  spherical  body,  we  have 


r(8)  -  sin  0, 


(3.2a) 


«(9)  -  1, 


(3.2b) 


s(0)  -  0 


(3.2c) 


This  also  includes  bodies  such  as  spherically  blunted  cones 
hemisphere-cylinder  bodies,  etc.  The  actual  requirement  here  is  that  the 
curvature,  K,  be  equal  to  unity  to  the  order  of  the  solution  obtained. 
The  results  of  the  second -truncation  problem  obtained  later  in  this 


chapter  are  applicable  to  any  body  for  which 
1  +  <2  sin  0  +  sin  0  +  ... 


0  where  x(s) 
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and,  of  course,  j  -  1  for  axisymmetric  flow. 

Since  there  Is  no  distinction  between  0  and  s  for  the  sphere, 
the  work  which  follows  will  be  given  in  terms  of  the  coordinate  s. 

The  shock  location,  A  ,  and  the  chock  angle,  0  ,  must  be  expanded 
in  order  to  be  applied  in  the  shock  boundary  conditions.  Thus 


A(s)  ■  sin  s  +  ... 

3 

0(s)  ■  0  -  (0X  8 in  s  +  03  sin  s  +  ...) 


(3.3a) 


(3.3b) 


If  0j^  *  03  ■  ...  «  0,  the  shock  will  be  a  sphere,  concentric  to  the 
body.  Thus  0^  ,  03  ,  , , ,  measure  the  deviation  of  the  shock  shape 
from  a  concentric  sphere. 

Equations  (3.1)  and  (3.2)  are  substituted  into  the  governing  partial 
differential  equations,  (2.2)-(2.7).  The  equations  are  expanded  into 
powers  of  sln(s),and  like  powers  of  sln(s)  are  equated.  The  lowest- 
ordered  terms  from  the  expanded  equations  yield 


2po 

v  p  +  p  v  -  77 —  (u,  -ev  ) 
oKo  Ko  o  1+en  1  0 

n  n 
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p  +  €  P  V  V  ■  0 
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n  n 
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(3.4b) 


(3.4c) 


H(T  )r  1  V>*  2 

V  p  -  p  v  T  -  - —  A  +  2e  T  +  - —  T  (3.4d) 

O  O  OOOCT  O  O  O  0 

n  n  _  nn  n  n 


p  -  p  T 
ro  Y  o  o 


(3.4e) 


Examining  these  equations,  we  find  that  they  involve  six  unknowns, 

u. ,  v  ,  p  ,  p»,  T  ,  and  p  ,  while  there  are  only  five  equations 
1  0  O  £  O  o 

available.  A  second  set  of  equations  can  be  obtained  by  equating  the 
next  higher -ordered  terms  of  each  expanded  equstlon.  This  results  in 
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(3.5b) 
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(3.5c) 
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There  are  now  a  total  of  ten  equations. but  the  number  of  unknowns  has 
become  eleven  as  u3  »  v2  *  ^2  ’  p2  *  ant*  ^4  have  been  added.  This 
pattern  continues  as  the  higher-order  equations  are  considered;  there 
is  always  at  least  one  more  unknown  than  there  are  equations.  Thus  we 
have  an  infinite  set  of  coupled  differential  equations.  In  the  case  of 
the  Blasius-series  solution  cited  earler,  the  coupling  between  the  equa¬ 
tions  goes  in  only  one  direction;  i.e.,  the  lower-order  terms  do  not 
depend  on  the  higher-order  terms,  and  thus  the  equations  can  be  solved  a 
few  at  s  time  until  as  many  terms  as  are  desired  have  been  computed.  In 
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Che  present  case,  however,  the  coupling  of  the  equations  Is  complete, 
and  it  would  appear  to  be  necessary  to  solve  all  of  the  equations 
simultaneously,  clearly  an  intractable  situation. 

By  introducing  approximations  for  the  excess  unknowns,  the  equa¬ 
tions  can  be  solved  to  give  approximate  values  for  the  coefficients. 

The  accuracy  of  these  solutions  will  depend  on  the  nature  of  the  approx¬ 
imations  Introduced  and  on  the  nature  of  the  problem  Itself.  However, 
the  sccuracy  of  these  approximations  can  be  systematically  evaluated 
and  improved. 

Consider  equations  (3.4)  again.  Only  one  of  the  second-order 
coefficients,  p2,  appears  and  it  occurs  in  only  one  term.  A  crude 
approximation  would  be  to  set  p2  =  0  .  This  would  permit  the  solution 
of  the  equations  and  yield  a  first  approximation  for  u^  ,  vq,  p^  , 
and  Tq  .  A  simultaneous  solution  of  equations  (3.4)  and  (3.5)  could 
si jo  be  obtained  if  p^  were  set  equal  to  zero.  This  second  approxima¬ 
tion  would  yield  not  only  an  approximate  solution  to  the  second-order 
coefficients  (u^  »  v2  ’  *  an<*  ^2^  ^Ut  wou^  also  yield  values 

for  the  first-order  coefficients  which  should  be  more  accurate  than 
those  previously  obtained.  The  Improved  accuracy  should  result  from 
the  fact  that  the  set  of  equations  (3.4)  would  be  solved  in  their  entirety 
with  more  accurate  values  of  p 2  than  were  used  in  the  first  approxi¬ 
mation.  This  would  then  afford  an  estimate  of  the  accuracy  of  the  first 
spproxlmation  (p2  ■  0)  .  This  process  could,  in  theory,  be  continued. 

At  each  step,  as  the  degree  of  approximation  is  improved  an  estimate 
of  the  accuracy  of  the  previous  approximation  is  obtslned. 
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The  process  of  setting  the  excess  unknowns  equal  to  zero  is  equiv¬ 
alent  to  assuming  that  the  variables  can  be  described  with  simple 
algebraic  expressions  rather  than  with  the  infinite  series  of  equations 
(3.1).  In  the  first  approximation  above,  a  one-term  expression  is  used, 
and,  in  the  second  case,  a  two -term  polynomial.  The  use  of  simple 
algebraic  expressions  to  evaluate  the  flow  around  a  blunt  body  has  been 
widely  employed  and  is  known  as  the  method  of  local  similarity  [19,  17, 
24,  28  ]•  The  series- truncation  method  may  be  considered  to  be  a 

generalization  of  the  method  of  local  similarity,  allowing  a  systematic 
improvement  of  the  degree  of  approximation.  As  such  the  series  trunca¬ 
tions  can  be  used  to  test  the  validity  of  the  concept  of  local  similar¬ 
ity. 

A  truncation  procedure  very  similar  to  the  scheme  outlined  above 
has  been  used  previously  [  31  ].  The  results  of  reference  31  indicate 
that  the  first  approximation  (p£  a  0)  of  this  scheme  would  yield  very 
poor  results.  In  order  to  improve  the  results,  it  would  then  be  neces¬ 
sary  to  retain  mor<:  terms  of  the  series  than  are  desirable  for  this 
problem.  Setting  p^  *  0  in  equation  (3.4)  is  equivalent  to  requiring 
a  zero  streamwise  pressure  gradient,  clearly  a  poor  approximation  for 
the  flow  around  a  blunt  body.  A  reasonable  value  for  the  pressure 
terms  is  needed,  and  this  can  be  accomplished  in  several  ways.  The 
"indented  truncation"  method  introduced  by  H.  C.  Kao  [21]  is  used  here. 
This  method  m-kes  use  of  equation  (3.5c)  to  approximate  the  value  of  p^. 

The  addition  of  equation  (3.5c)  to  equations  (3.4)  results  in  eight 
unknowns  appearing  in  six  equations.  The  "excess"  unknowns  which  are 
truncated  to  permit  a  solution  are  P2  and  V2*  The  solution  obtained 


-36- 


Therefore  this 


in  this  manner  provides  a  reasonable  value  for  p2 
first  approximation  is  far  more  accurate  than  the  one  described  earlier. 
Equations  (3.4)  and  equation  (3.5c)  (with  p2  -  v2  -  0)  constitute  the 
first- truncation  problem.  The  values  of  the  variables  which  are  obtained 
from  the  first  truncation  are  denoted  by  v^  and 

(1)  1  o  ’  °  *  ° 

P2  .  The  boundary  conditions  for  the  first  truncation  are  obtained 
from  the  expansion  of  equations  (2. 10) -(2. 17)  by  use  of  equations  (3.1)- 
(3.3).  The  boundary  conditions  which  apply  to  the  variables  of  the 
first  truncation  are  given  below.  At  the  body  surface  the  velocity- 
and  temperature-jump  conditions  yield 


M.CT  CO)) 
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(3.6a) 


°  Po(0) 


v? 


ec  li(T  (0) ) 

T-(°)  *  + - /  - 1^(0)  T^  (0)  (3.6b) 


o'  o 


vq(0)  =•  0 


(3.6c) 


where  the  specified  body  temperature  has  been  expanded  into 


Tb(s) 


+  T,  sin  s  + 


(3.7) 


Immediately  behind  the  shock, the  following  conditions  must  be  applied. 
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"i<v  ■ 1  ■  v^i  [i  ■  i2]”1 


(3.8a) 


e  v  (A  )  -  ^  + - - — r 

Y+l  (Y+1)M? 


(3.8b) 


p  (A  )  *  - 

°  0  Y+l 


L  2V  J 


(3.8c) 
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(3.8d) 


(3.8e) 


The  following  geometric  relation  between  0^  and  ^  is  obtained 
by  substituting  equations  (3.3)  into  equation  (2.19). 


2  e 


1  l+e  A 


(3.9) 


The  solution  to  this  first  truncation  problem  is  discussed  in  sec¬ 
tion  C  where  it  will  be  found  that  it  is  also  necessary  to  approximate 
the  value  of  0^  in  equations  (3.8a)  and  (3.8d). 

To  assess  the  accuracy  of  the  first  truncation  (and  to  simultan¬ 
eously  Improve  it),  equations  (3.4)  and  (3.5)  are  solved  simultaneously 
by  making  a  suitable  assumption  for  the  value  of  p^  .  As  was  the  case 
in  the  first  truncation,  the  necessary  value  is  obtained  from  the  normal 
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momentum  equation.  The  third-order  terms  of  the  expanded  normal  momen¬ 
tum  equation,  when  equated,  yield 
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As  before,  it  Is  necessary  to  truncate  two  excess  unknowns.  Therefore 

we  equate  p^  and  to  zero.  Equations  (3.4),  (3.5)  and  (3.10) 

(with  P4  "  =  0  )  constitute  the  second -truncation  problem, and  the 

(2)  (2) 

values  obtained  from  Its  solution  are  denoted  by  u'  ,  v  , 

1  o 

(2)  (2)  (2)  (2) 

u^  ,  v2  ,  ^2  *  ^4  *  *n  theory*  this  procedure  could  be  continued 

with  more  equations  being  solved  at  each  step.  As  the  order  of  the 
truncation  is  Increased  in  this  manner,  the  problem  more  closely 
approximates  the  infinite  set  of  coupled  differential  equations  which 
govern  the  flow  variables.  It  is  assumed  that  the  solutions  to  the 
truncated  equations  approach  the  solution  of  the  full  equations.  That 
is,  the  basic  premise  of  the  method  of  series  truncation  is  that  the 
sequence  of  solutions  for  each  variable,  e.g.  uj^ ,  u|^ ,  ...  u|n^ ,  ... 
converges  to  the  desired  solution  of  the  infinite  set  of  ordinary 
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differential  equations.  A  brief  examination  of  the  increasing  complexity 
of  equations  (3.4),  (3.5)  and  (3.10)  show«  that  if  the  method  is  to  be 
practical  the  convergence  of  the  sequence  of  solutions  must  be  rather 
rapid.  To  proceed  beyond  the  second  truncation  becomes  rather  difficult 
for  this  problem.  However,  there  is  reason  to  believe  that  the  conver¬ 
gence  is  rapid  for  this  particular  problem.  In  fact, the  method  of  local 
similarity  assumes  that  the  flow  is  adequately  described  by  expressions 
equivalent  to  those  used  in  the  first  truncation. 

To  complete  the  formulation  of  the  second  truncation  problem, 
additional  boundary  conditions  are  obtained  in  the  same  manner  that 
equations  (3.6),  (3.8)  and  (3.9)  were  obtained.  At  the  wall,  we  have 


v2(0)  -  0 


(3.11c) 
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After  considerable  algebraic  manipulation  the  following  boundary  con¬ 


ditions  at  the  shock  are  obtained. 
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These  boundary  conditions  are  used  together  with  those  previously  given 
in  equations  (3.6)-(j.9)  to  complete  the  formulation  of  the  second 
trur.ction  problem. 


e[v2<Ao)  +  *2  VV]  -  - 


/  1  N 

^  '  M2  > 


(3.12b) 


-41- 


To  evaluate  the  shear  stress  and  heat  transfer  at  the  body  surface 
from  these  series  solutions,  It  Is  necessary  to  expand  the  expressions 
for  shear  stress,  t  ,  and  heat  transfer  rate,  q  ,  by  substituting 
equations  (3.1)  Into  equations  (2.22).  The  shear  stress  Is  then  given 
by 

t(b,0)  ■  e(T^  sins  +  sln^s  +  ...)  (3.14) 

where 


for  Tl  ,  t3 


f 


The  solution  of  the  first  and  second  truncations  will  He  obtained 
for  several  blunt-body  flows  In  section  D.  These  solutions  provide  a 
measure  of  the  range  of  validity  of  the  flow  model  adopted  In  Chapter 
II  and  provide  a  test  of  the  validity  of  the  concept  of  local  similarity 
as  well  as  fulfilling  their  original  purpose  of  providing  Initial  data 
at  the  axis  of  symmetry.  Before  the  equations  are  actually  solved, 
however,  some  of  the  previous  Investigations  based  on  series  truncations 
will  be  discussed  In  section  B. 

B.  Previous  Applications  of  Series  Truncations 

There  Is  no  mathematical  proof  of  the  convergence  of  the  method 
described  In  the  previous  section  and  thus  no  assurance  of  Its  validity. 
Further,  the  necessity  of  obtaining  a  rapid  convergence  means  that  an 
evaluation  of  the  effectiveness  of  the  method  must  come  from  an  analysis 
of  the  application  of  the  method  to  specific  examples.  Most  of  the  pre¬ 
vious  applications  of  series  truncations  have  been  to  the  blunt-body 
problem  [31,  6,  7,  21,  36,4].  However,  other  examples  of  its  use  include 
Incompressible  boundary- layer  flow  over  paraboloids  [  10  ];  viscous  flow 

4 

over  a  semi-infinite  plate  (  8  ];  and  Incompressible,  viscous  flow  around 
a  circular  cylinder  [34  ].  These  Investigations,  except  for  ref.  10,  have 
been  reviewed  by  M.  Van  Dyke  In  reference  35.  The  work  that  Is  discussed 
in  this  reference  clearly  illustrates  certain  Important  aspects  of  the 
series- truncation  method,  In  particular  It  should  be  noted  that  the 
accuracy  of  the  method  may  be  quite  dependent  upon  the  form  of  the  expan¬ 
sion  and  that  the  method  may  be  expected  to  work  best  for  problems  in 
which  the  flow  field  Is  thin  (one  dimension  much  smaller  than  the 

other).  This  later  aspect  Is  anticipated  In  the  use  of  the  method  of 
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local  similarity  for  such  problems. 

The  coupling  of  the  lowest-ordered  terms  to  the  higher-ordered 
terms  is  a  consequence  of  the  upstream  influence  which  occurs  in  the 
flow.  That  is,  it  represents  the  fact  that  the  flow  at  the  axis  depends 
upon  the  flow  in  the  region  away  from  the  axis.  Clearly  the  method  will 
be  most  successful  when  this  dependence  is  weak.  This  will  be  the  case 
when  the  equations  are  "nearly  parabolic"  (or  hyperbolic)  and/or  the 
region  of  Interest  is  thin.  For  thin  regions  the  solution  should  depend 
most  strongly  on  the  locally  imposed  boundary  conditions.  In  fact,  in 
the  limit  as  the  width  of  the  region  goes  to  zero,  e.g. ,  in  boundary- 
layer  theory  and  in  thin-shock-layer  theory,  the  governing  differential 
equations  become  parabolic.  There  is  then  no  upstream  influence  at  all, 
and  no  truncations  are  necessary. 

It  can  be  seen  that  the  investigations  cited  above  pose  tests  of 
varying  degrees  of  difficulty.  The  boundary- layer  flow  over  a  para¬ 
boloid  considered  by  R.  T.  Davis  in  reference  10  is  exactly  the  type  of 
flow  for  which  the  method  is  expected  to  work.  Since  the  equations  are 
parsbolic,  an  expansion  about  the  axis  of  symmetry  would  need  no  trun¬ 
cation  snd  would  lead  to  the  Blaslus-serles  solution.  However,  Davis 
has  used  a  variation  of  the  method  which  is  referred  to  as  the  method  of 
local  truncations.  The  variables  are  expanded  about  sn  arbitrary  down¬ 
stream  station.  The  resulting  problem  is  then  solved  st  various  down¬ 
stream  positions,  and  the  results  of  each  solution  sre  used  only  locally. 
A  more  accurste  representation  of  the  downstream  flow  can  be  obtained  in 
this  manner  than  could  be  obtained  from  s  reasonable  number  of  terms  of 
an  expsnsion  at  the  sxis.  However,  this  method  has  two  disadvantages: 
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the  flow  symmetry  is  lost,  thereby  requiring  more  terms  in  the  expansion; 
and  the  coupling  between  different-ordered  terms  must  now  describe  the 
influence  of  the  flow  in  regions  both  upstream  and  downstream  of  the 
point  of  expansion.  In  the  flow  considered  by  Davis,  there  is  no 
upstream  inf luence, but  a  coupling  does  occur  since  the  solution  at  the 
expansion  station  depends  on  the  flow  near  the  axis.  However,  for  the 
bound ary -layer  equations,  it  is  known  that  this  dependence  is  weak  since 
the  Influence  of  the  initial  data  dies  out  rapidly  away  from  the  initial 
station.  Even  though  this  flow  problem  is  a  weak  test  of  the  method,  it 
is  a  significant  test  since  exact  numerical  solutions  are  available  from 
finite-difference  methods.  Davis  shows  that  the  second  truncation  pro¬ 
duces  highly  accurate  results  over  the  entire  paraboloid. 

The  investigation  of  reference  8  uses  local  truncations  to  solve 
the  Navier-Stokes  equations  for  flow  over  a  semi-infinite  plate.  Since 
the  equations  are  elliptic,  this  problem  poses  a  more  difficult  test  of 
the  method.  Although  no  exact  solution  exists  for  comparison,  the 
results  obtained  for  this  problem  are  remarkable  and  substantiate  the 
value  of  the  method. 

Van  Dyke  has  posed  a  particularly  severe  test  of  the  method  of 
series  truncation  in  reference  34  .  He  has  used  both  the  linearized 

Oseen  equations  and  the  full  Navier-Stokes  equations  to  describe  the 
flow  around  a  circular  cylinder.  The  method  is  not  expected  to  be 
highly  accurate  for  such  elliptic  equations  in  a  thick  (infinite)  region 
and.  Indeed,  the  second-truncation  solution  of  the  Oseen  equations  shows 
only  qualitative  agreement  with  the  exact  Oseen  solution.  However,  Van 
Dyke  has  shown  that  a  knowledge  of  this  error  when  combined  with  the 
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first-truncation  solution  to  the  Navler-Stokes  equations  yields  very 
useful  results  even  for  this  particularly  severe  problem. 

The  analyses  of  the  blunt-body  problem  based  on  series  truncations 
are  also  encouraging.  Swigart  [  31]  solved  the  invlscid  flow  of  a  per¬ 
fect  gas  over  spheres  and  paraboloids  at  small  angles -of-attack.  Three 
and  four  truncations  produced  accurate  results  at  the  axis.  However,  at 
the  sonic  line  the  results  appeared  to  be  considerably  less  accurate. 

R.  J.  Conti  [6,7],  investigating  the  invlscid  flow  of  a  chemically 
reacting  gas,  and  H.  C.  Kao  [21],  investigating  the  viscous  shock  layer, 
improved  the  convergence  of  the  method  so  that  reasonable  results  were 
obtained  near  the  axis  with  only  two  truncations.  This  is  particularly 
Important  for  the  viscous  shock  layer  since  it  is  difficult  to  proceed 
beyond  the  second  truncation  due  to  the  increasing  complexity  of  the 
equations.  Van  Dyke  [  36  ]  re-examined  Swigart' s  solution  for  a  para¬ 
boloidal  shock  wave  and  showed  that  a  modified  expansion  procedure  would 
produce  a  solution  of  remarkable  accuracy  over  the  entire  body.  P.  Cheng 
and  W.  Vincentl  [  4  ]  have  used  this  knowledge  to  Investigate  the  flow  of 
a  radiating,  invlscid  gas  over  a  paraboloid.  Each  of  the  above  analyses 
was  based  on  the  inverse  method  in  which  the  shock  shape  is  given  and 
the  body  which  produced  that  shock  shape  is  calculated.  This  simplifies 
the  computational  procedures  but  leaves  some  ambiguity  as  to  what  prob- 
blem  has  been  solved,  especially  for  the  flow  of  a  viscous,  heat  conduct¬ 
ing  fluid.  The  investigation  by  Kao  is  of  particular  interest  for  the 
work  that  follows  in  this  chapter. 

The  success  of  the  method  of  series  truncation  in  the  analysis  of 
the  difficult  problems  described  above  indicates  that  it  is  a  valuable 


analytical  method  and  justifies  its  use  for  the  computation  of  the 
flow  variables  at  the  axis  of  symmetry  of  a  blunt  body. 

C.  Solution  of  the  Equations 

The  set  of  equations,  (3.4)  and  (3.5c),  which  constitutes  the  first- 
truncation  problem, is  a  seventh-order  system  of  nonlinear  ordinary  dif¬ 
ferential  equations  for  the  six  unknown  functions,  u^,  vq,  pq,  Tq 
and  Equation  (3.4e)  may  be  used  to  eliminate  pQ  immediately. 

Together  with  the  boundary  conditions  (equations  (3.6),  (3.8)  and  (3.9)), 
this  system  of  equations  comprises  a  two-point  boundary-value  problem. 

Due  to  the  cortplexity  of  these  equations,  solutions  must  be  obtained  by 
numerical  methods.  The  equations  are  solved  numerically  by  reducing  the 
two-point  boundary-value  problem  to  an  equivalent  initial-value  problem. 
This  is  accomplished  by  guessing  the  values  of  the  unknown  variables  at 
one  of  the  boundaries.  The  equations  may  then  be  integrated  to  the  other 
boundary  by  using  any  of  the  standsrd  integration  procedures  designed 
for  use  with  automatic  digital  computers.  The  computations  of  this 
chapter  were  made  on  the  Burroughs  B5500  computer  using  the  Kutta-Merson 
method.  The  guessed  initial  values  must  then  be  altered  until  the 
boundary  conditions  at  the  second  boundary  are  satisfied.  This  iteration 
on  the  guessed  initial  values  is  accomplished  by  using  the  Newton-Raphson 

l 

method  to  calculate  the  necessary  corrections  to  the  initial  values. 

If  the  original  guessed  values  are  not  sufficiently  close  to  the  correct 
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A  description  of  this  method  may  be  obtained  from  any  standard  text 
on  numerical  methods,  e.g.,  Introduction  to  Numerical  Analysis,  F.  B. 
Hildebrand,  McGraw-Hill  Book  Company,  Inc.,  1956. 


values,  the  Iterations  based  on  the  Newton-Raphson  method  may  not  con¬ 
verge  or  may  converge  rather  slowly.  However,  in  such  cases,  it  has 
been  found  that  convergence  can  be  obtained  by  using  only  a  fraction  o i 
the  corrections  predictedby  the  Newton-Raphson  method.  The  exact  value 
of  this  fraction  is  selected  so  that  the  errors  in  the  boundary  condi¬ 
tions  are  minimized  at  the  second  boundary.  The  remainder  of  this 
section  is  a  description  of  the  details  of  these  computations.  The 
reader  is  referred  to  section  D  for  the  results  of  these  computations. 

In  Chapter  II  it  was  noted  that  difficulties  with  this  computational 
method  may  be  expected  when  the  body  surface  is  a  characteristic  sur¬ 
face  for  the  governing  partial  differential  equations.  These  difficul¬ 
ties  are  also  Illustrated  in  the  remainder  of  this  section. 

The  application  of  any  of  the  numerical  techniques  for  Integration 
requires  the  reduction  of  the  set  of  equations  to  a  system  of  first-order 
differential  equations  having  the  form 

^  -  f(n,y)  (3.16) 

where  y  and  f  are  vector  quantities.  This  reduction  i3  easily 

accomplished  by  considering  u  and  T  to  be  separate  unknowns  and 

in  °n 

d(ux)  d(To) 

introducing  two  additional  equations,  -jj —  ■  u.^  and  —  3  Tq  . 

n  n 

The  vector  unknown,  y,  has  seven  cosponents  for  this  first -truncation 

problem:  y  -  (u^  ,  vq,  pQ,  p2,  Tq,  Tq  ).  In  order  to  obtain  the 

n  n 
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dp  dv 

first  derivatives  •“  and  -r-2-  explicitly,  as  required  in  (3.16), 

dn  dn 

it  is  necessary  to  algebraically  combine  the  continuity  (3.4a),  normal 
momentum  (3.4c)  and  state  (3.4e)  equations.  The  remaining  first  deriv- 
tlves  are  obtained  immediately  from  the  appropriate  equations.  The 
seven  components  of  equation  (3.16)  for  the  first- truncation  problem 

l 

are 


dul 

~  fl(n.y)  -  U1 


(3.17a) 


du 


IT  ■  ¥"->•>  ■  J<t7 


(A..  £V< 

l  1+en  V°U^-n  1+en 


2Ji L 

1+en 


dv 


n'<V 

2  £ui  ■  JxtT  <U1  ■  '"P  To 

n  **  o  n  r 


(3.17b) 


v  T 

7  0  0 
1+fn  <U1  ’  €Vo>  +  -T-11 


d r  -  f3<n’*>  - 


(3.17c) 


1  -  € 


Y-l  T  '■ 
'  o 


dPn  2 

dir  “  f4(n’y)  =  £  poVo  f3(n*y) 


(3 . 17d) 


dP? 

d r  -  f5<n*y>  - 
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eu.  v 
1  o 


2  n 


I+ciT"  "  e  vo  f3(n->r) 


1+en 


J  (3. 17e) 


The  superscript  (1)>  denoting  the  values  obtained  from  the  first 
truncation,  has  been  dropped  for  simplicity. 
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where 


dT 

dT  ■  £6(n'y)  *  Tc 


(3. 17f> 


ST1  -  v„  ’  *>o  Tol  ~ 

-  T  <V  1 

T  2e  + - 2-  T 

°“  L  H(TJ  °nj 


(3.17g) 


Y  P£ 

-v  m  •—  — ■ 

o  Y“1  T 


There  are  nine  boundary  conditions  available  in  equations  (3.6), 
(3.8)  and  (3.9).  However,  there  are  three  additional  unknown  quantities, 
Aq  ,  and  0^  ,  which  appear  in  these  boundary  conditions.  Thus 

there  is  actually  one  less  boundary  condition  available  than  is  needed  to 
define  a  solution  to  equations  (3.17).  The  problem  is  made  determinate 
by  truncating  the  series  that  describes  the  shock  position  or  shock 
angle,  equation  (3.3).  Setting  0^  ■  0  (and  therefore  Aj  “  0  )  is 
equivalent  to  assuming  that  the  shock  and  body  are  concentric.  Thue, 
with  this  approximation,  the  first  truncation  problem  becomes  consistent 
with  the  analyses  based  on  local  similarity  and  the  thin-shock-layer 
model. 

Equations  (3.6  a,b,c)  provide  three  relations  for  the  seven  depen¬ 
dent  variables  at  the  body  surface.  Thus, to  begin  the  integration  of 
equations  (3.17)  at  the  wall, it  is  necessary  to  guess  the  values  of  four 
of  the  variables.  The  recommended  procedure  is  to  choose  values  of 
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P_(0)  *  P,(0)  I  u-  (0)  aM  T  (0).  Equation  (3.6c)  gives  v  (0)  -  0, 

o  Z  I  o  0 

n  n 

and  equations  (3.6a)  and  (3.6b)  must  be  solved  to  obtain  the  values  of 
u^(0)  and  Tq(0)  that  are  consistent  with  the  guessed  values  for 
P  (0),  u.  (0)  and  T  (0).  The  solution  of  equation  (3.6b)  for  T  (0) 

O  1  O  0 

n  n 

Is  not  Immediate  since  Tq  appears  In  an  Implicit  manner.  Using  a 
power  law  for  the  viscosity , 


TU)  , 


we  may  rewrite  equation  (3.6b)  In  the  form 


v°> 


(3.18) 


To  <°> 

1  -  ec 

1  Y  PQ(0) 


To(°) 


0)  -  j 


This  relation  may  be  solved  by  successive  substitution.  That  Is,  sub¬ 
stituting  an  aoproximate  value  of  T  (0)  ,  say  T  (0)  «  T.  ,  into 

O  0  t> 

0 

the  right-hand  side  yields  an  Improved  value  for  Tq(0)  .  This  value  is 

in  turn  substituted  into  the  right-hand  side.  This  process  is  continued 

until  convergence  Is  obtained.  Note  that  for  the  special  esse,  u>  *  1/2, 

the  solution  is  obtained  Immediately  since  T  (0)  is  removed  from  the 

0 

right-hand  side  of  the  equation.  With  Tq(0)  available,  equation  (3.6a) 
can  be  solved  for  u^(0).  The  functions  f(n,y),  defined  by  equations 
(3.17),  may  now  be  evaluated  at  n  ■  0,  and  the  numerical  Integration 
of  equations  (3.17)  can  be  performed  In  a  step-wise  manner  until  the 
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shock  is  located.  The  shock  condition  on  v  has  been  found  to  work 

o 

well  for  the  determination  of  Aq.  The  numerical  integration  is  con¬ 
tinued  for  increasing  n  until  the  value  of  vq  satisfies  equation 

(3.8b).  This  defines  the  value  of  A  .  There  remain  four  shock  con- 

o 

ditions  which  are  satisfied  by  Iterating  on  the  four  initial  values 

u.  (0)  ,  T  (0)  ,  p  (0)  and  p  (0)  .  Due  to  the  nonlinearity  of  the 
n  n  z 

equations,  several  such  iterations  using  the  Newton-Raphson  method  are 
generally  required  to  obtain  the  correct  values  of  the  variables  at  the 
wall.  Each  Iteration  proceeds  in  the  following  manner. 

1)  Each  of  the  guessed  initial  values  is  varied  by  a  small 
Increment.  The  equations  are  integrated, and  the  effect 
of  the  change  in  the  initial  value  on  the  unsatisfied 
shock  conditions  is  evaluated.  This  step  requires  the 
numerical  integration  of  equations  (3.17)  four  additional 
times. 

2)  With  the  assumption  that  the  shock  conditions  depend  on 
changes  in  the  initial  data  in  a  linear  manner,  a  correction 
to  the  initial  data  is  computed. 

3)  The  equations  are  then  integrated  with  the  corrected 
initial  values, and  the  shock  conditions  are  evaluated  again. 

This  process  continues  until  the  shock  conditions  are  satisfied  to  the 
desired  degree  of  accuracy. 

The  application  of  the  temperature- jump  condition,  equation  (3.6b) 

or  (3.18),  is  simpler  if  an  alternative  procedure  of  guessing  Tq(0) 

and  computing  the  consistent  value  of  T  (0)  is  adopted.  This 

n 


procedure  works  quite  well  for  low  values  of  the  Reynolds  number. 

However,  numerical  difficulties  frequently  occur  when  e  Is  very  small 

(large  Reynolds  number).  The  computed  value  of  T  (0)  (and  therefore 

o 

n 

the  solution  to  equations  (3.17))  Is, In  this  case, quite  sensitive  to 

the  guessed  value  of  Tq(0).  Hence,  In  the  Iteration  procedure 

described  above,  the  Increments  added  to  T^(0)  must  be  kept  very  small. 

It  occasionally  becomes  difficult  to  keep  the  errors  that  result  from 

the  numerical  Integration  smaller  than  this  Increment  in  Tq(0).  If 

this  is  not  accomplished,  the  predicted  corrections  to  the  Initial 

conditions  are  meaningless,  and  the  iterations  do  not  converge.  For  this 

reason  It  was  recommended  that  the  value  of  T  (0)  be  guessed. 

n 

Probstsln  and  Kemp  [  28  ]  considered  several  approximate  methods 

based  on  local  similarity  and  showed  that  some  of  the  schemes  could 

lead  to  an  overdetermined  system  of  differential  equations.  Howe.er, 

one  generally  encounters  a  system  which  appears  to  be  underdetermined. 

In  the  case  of  the  first* truncation  problem  described  above,  this  was 

resolved  by  requiring  the  body  and  shock  to  be  concentric.  In  the  thin- 

shock-layer  investigations  based  on  local  similarity,  e.g.  ref.  19  , 

the  systems  of  ordinary  differential  equations  usually  appear  to  be 

underdetermined  even  though  such  investigations  assume  that  the  body  and 

shock  are  concentric.  This  is  a  result  cf  the  fact  that  the  (uv  V 

n  n 

term  has  been  retained  in  the  equations.  The  order  of  the  system  Is 
then  one  greater  than  the  system  defined  by  equations  (3.17).  Hence, 
one  additional  boundary  condition  or  constraint  is  required  to  determine 
the  value  of  v  (0).  It  was  also  noted  In  Chapter  II  that  keeping 
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(uv  )  In  the  normal  momentum  equation  could  lead  to  difficulties  in 
n  n 

the  formulation  of  the  problem.  Thus  it  is  of  interest  to  examine  the 

form  that  equations  (3.17)  would  take  if  the  (uv  )  terra  were  included. 

~  n  n 

In  this  case,  the  vector  unknown,  y  ,  would  have  eight  components 

since  v  would  have  to  be  considered  as  a  separate  variable.  Several 
n 

of  the  equations  in  (3.17)  would  change  since  the  normal  momentum  equa¬ 
tion  is  of  the  form 


2 

e  P, 


v  v 
o  o  o 


-  *<V 

nn 


dv  dp 

In  particular,  the  expressions  for  and  ^2  given  in  (3.17c)  and 

(3.17d)  would  be  replaced  with 


dv 
_ o 

dn 


=  f3(n,y) 


\ 

k 


.i 


and 


dp 

if  '  £4<n>1,) 


p  X  p 

ro  o  ro 

— - — -  +  — 

T  v 

o  o 


2(u  -  ev  ) 

1  o 

1+tn 


The  eighth  equation  would  be 


dv 

-if  ■  ■  •  jhj 


V  V 

0  o 


1 

nJ 


1 

The  right-hand  side  of  this  equation  is  not  complete,  but  in  this 
form  the  equation  has  the  same  essential  features  as  the  more  complete 
eouation. 
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This  is  consistent  with  the  analysis  in  Chapter  II  based  on  the  char¬ 
acteristic  surfaces  of  the  governing  partial  differential  equations 

(page  22  ).  There  it  was  noted  that  the  inclusion  of  the  (uv  )  tern 

r  n  n 

meant  that  streamlines  were  characteristics.  Thus,  for  no  mass  transfer, 
the  body  surface  is  characteristic,  and  one  cannot  determine  the  highest- 
ordered  derivatives  normal  to  the  body  unless  the  specified  initial 
data  satisfy  a  compatabllity  equation.  It  is  now  seen  that  the  compat¬ 
ibility  equation  requires 


2u.(0)  -  v  (0)  ■  0  . 
i  o 

n 

In  the  case  of  no  slip  at  the  boundary,  this  simplifies  to  a  requirement 

that  v  (0)  *  0.  This  compatabllity  condition  is  the  additional  con- 
n 

stralnt  or  boundary  condition  which  is  used  in  numerous  investigations 
[  19,  21,  24,  28  ]  to  make  the  system  of  equations  determinate. 

However,  it  was  noted  in  Chapter  II  that  all  difficulties  are  not 
removed  by  this  step  as  there  is  no  guarantee  that  the  solution  is 
unique  when  the  initial  data  satisfy  this  compatabllity  relation.  Shlh 
and  Krupp  [30  ]  have  computed  an  alternative  solution  to  one  of  the 

-55- 


g 


!  ' 


j 


_ 


f 

ti 


examples  considered  by  Ho  and  Probstein  [  19 J,  thus  verifying  the  non¬ 
uniqueness  of  such  a  formulation  of  the  problem. 

In  an  investigation  similar  to  that  of  Ho  and  Probstein,  L.  Goldberg 
[  17]  has  considered  examples  in  which  there  is  mass  transfer  at  the 
wall.  In  such  a  case  there  is  no  comparability  relation  at  the  wall,  and 
the  system  of  differential  equations  again  appeara  to  be  underdetermined. 
Goldberg  uses  the  Integral  form  of  the  mass -conservation  equation  (c.f. 
equation  (2.18))  as  the  necessary  additional  constraint.  However,  in 
general,  this  is  not  an  independent  condition.  If  the  continuity  equa¬ 
tion  is  satisfied  throughout  the  shock  layer  and  all  the  boundary 
conditions  are  properly  enforced,  the  overall  mass  conservation  is 
assured.  An  examination  of  Goldberg's  equations  shows  that  in  the  con¬ 
tinuity  equation  a  factor  of  the  form  1+en  has  been  approximated  by 
unity.  This  is  consistent  with  the  thin-shock-layer  model,  but  a  simi¬ 
lar  approximation  in  the  integral  equation  has  not  been  made.  Thus  the 
relation  is  independent  due  to  the  approximate  nature  of  the  continuity 
equation  that  is  used.  The  problem  is  then  determinate,  and  the  shock 
stand-off  distance  can  be  computed  (a  concentric  shock  interface  and 
body  surface  had  already  been  assumed). 

As  noted  in  Chapter  II,  Shih  and  Krupp  modify  their  governing 

equations  in  such  a  way  that  the  special  role  of  the  (‘jivn)n  term  is 

altered.  In  this  manner  the  difficulties  associated  with  obtaining 

v  (0)  are  apparently  avoided  (although  this  point  is  not  clear).  How- 
o 

n 

ever,  they  have  not  assumed  a  concentric  shock  and  body,  thus  raising  by 
one  the  number  of  boundary  conditions  needed.  This  additional  relation 
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Is  again  taken  to  be  the  integral  form  of  mass  conservation.  It  is 
not  clear  why  this  equation  should  give  an  independent  constraint  for 
their  analysis,  and  this  point  is  not  discussed  by  the  authors.  However, 
it  is  probably  related  to  the  fact  that  the  shock  conditions  are  applied 
only  at  a  specified  number  of  stations  downstream  of  the  axis  of  symmetry. 

It  is  quite  obvious  that  in  any  analysis  of  the  type  considered  in 
this  chapter  the  question  of  the  determinacy  of  the  system  of  differ¬ 
ential  equations  may  be  rather  complicated.  For  the  present  analysis  it 
appears  to  be  necessary  to  make  some  assumption  about  the  shock  slope. 
Hence  0^  was  taken  to  be  zero  in  the  first -truncation  problem 
described  above.  The  use  of  the  Integral  form  of  mass  conservation 
gains  nothing.  If  equation  (2.18)  is  expanded,  it  becomes,  for  the  first- 
truncation  problem  for  flow  around  a  sphere, 

^o 

(l+eAQ) 2  =  2  J  p^d+en)  dn  .  (3.19) 

o 

If  the  continuity  equation  (3.4a)  is  rewritten  in  the  form 

((1+en)2  povo)n  -  p^U+en) 

and  substituted  into  equation  (3.19),  equation  (3.19)  immediately  reduces 
to  po(Ao)vo(Aq)  -  1.  This  result,  however,  has  already  been  assured  by 
the  application  of  the  Ranklne-Hugoniot  shock  conditions. 

The  computation  of  the  second -truncation  problem  proceeds  in  the 
same  manner  as  did  the  first- crunc at ion  problem.  The  vector  unknown, 
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y  ,  of  equation  (3.16)  now  has  thirteen  cor^onents:  y  ■  (u  ,  u  ,  v  , 

1  ln  0 


po’ 


P2’  Td’  Tc 


3,  u3  ,  v2>  p^,  T2,  T2  ).  Equations  (3.4),  (3.5)  and 
n  n 


(3.10)  are  reduced  to  a  system  of  thirteen  first-order  differential 

equations.  Of  the  first  seven  of  these  equations,  all  except 
dp2 

*  f^(n,y)  are  identical  to  the  corresponding  equations  of  the  first 

dp2 

truncation  problem.  must  now  be  obtained  from  an  algebraic  com¬ 

bination  of  equations  (3.5a),  (3.5c)  and  (3.5e).  Thus  the  explicit 
formulation  of  the  second- truncation  problem  consists  of  equations 

l 

(3. 17a) -(3. 17d) ,  (3.17f),  (3.17g)  plus  the  following  equations. 


2  2 

dP?  (  6  \  9 

d/  =  f5(n,y)  =  1  +| - 222— J  ^(n)  "  6  v^n)  (3<20) 


-  6  P  V 
O  ^O  01 


where 


«oul 


Pi*0'  =  ufe  [“l  +  s(2v2  •  VD  -  6  [po(v2  -  V  +  p2Vo]  £3<"’y) 


and 


‘  ~  i (  P2  To  ‘  T2  V“'y) 

Y-l  T  \  n  n/ 


’  2v  T 
o 


"T - f3(n'y)  *  P0V2  p" 


L°o(2u3  ‘  ev2)  +  c2(2u1  ‘  fV] 


The  superscript  (2),  denoting  the  values  obtained  from  the  second 
truncation, have  been  dropped  for  simplicity.  It  should  be  remembered, 
however,  that  these  variables  are  a  second  approximation  t  >  the  correct 
values  of  the  flow  variables.  In  particular  it  should  be  noted  that  the 
function  f.(n,y)  which  appears  in  equations  (3.21c)  and  (3.21f)  is 
3  (2) 

the  function  f'  7  defined  in  equation  (3.20)  and  not  the  function 
(1  \  5 

f^  '  defined  previously  in  equation  (3.17e). 
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Y  p0  Y 

P  — r  z~  and  p-  ■  — - 

Ko  y-1  Tq  k2  Y“1 


p2  To  -  Po  T2 


Despite  the  complexity  of  these  expressions,  they  are  in  a  form 
which  is  easily  evaluated  numerically  on  a  computer,  given  a  value  of 
n  and  the  corresponding  values  of  the  dependent  variables.  There 
are  seventeen  boundary  conditions  available  from  equations  (3.6),  (3.8), 
(3.9),  (3. 11) -(3. 13).  The  conditions  at  the  shock  contain  an  additional 
five  unknowns,  Aq,  Aj,  A^,  0 1  and  03>  which  describe  the  shock 
shape.  Thus  the  problem  is  again  underdetermined,  and  some  suitable 
assumption  regarding  the  shock  shape  must  be  made.  Unlike  the  first- 
truncation  problem,  setting  A^  ■  0  does  not  imply  that  03  ■  0. 

Hence  we  can  truncate  either  the  shock-position  series  or  the  shock- 
angle  series  with  equation  (3.13)  supplying  the  coefficient  of  the 
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other  series.  Since  there  Is  no  clear  choice  here,  It  also  seems 
reasonable  to  truncate  both  series,  In  which  case  the  geometric  relation 
given  by  equation  (3.13)  is  not  satisfied.  These  various  possibilities 
are  considered  in  the  next  section.  It  should  be  not'd  that  0^  and 
aTe  computed  in  this  second -truncation  problem.  Therefore  a  check 
is  provided  on  the  accuracy  of  the  assumption  that  0^  ■  “  0  in 

the  first- truncation  problem. 

Now  it  is  necessary  to  guess  seven  unknown  values  at  the  body 

surface:  Uj^  (0),  pQ(0),  p2(0),  Tq  (0),  u3(0),  p^(0)  and  T2  (0).  The 
n  n  n 

boundary  conditions  (3.6a-c)  and  (3.11a-c)  provide  the  remaining  six 

values.  The  equations  are  then  integrated  numerically  until  the  value 
of  vq  satisfies  equation  (3.8b).  This  defines  the  value  of  Aq.  The 
value  of  0^  can  be  obtained  by  solving  equation  (3.12b).  Since  this 
equation  is  quadratic  in  0^  ,  two  values  are  obtained.  However,  it 
har  been  found  in  all  of  the  examples  considered  thus  far  that  the 
larger  value  of  0^  corresponds  tc  a  shock  of  negative  curvature 
(  ^  <  0)  and  thus  is  unacceptable.  The  requirement  for  a  positive 
shock  curvature  is  found  by  differentiating  equation  (3.3b).  At  the 
axis,  a  positive  curvature  requires  that  <  1  .  After  Aq  and  0^^ 
are  determined,  there  are  seven  boundary  condltlons--equations  (3.8a), 
(3.8c)-(3.8e) ,  (3.12a),  (3.12c),  (3. 12d) — which  can  be  used  to  determine 
the  correct  values  of  the  seven  guessed  initial  conditions.  The  itera¬ 
tion  procedure  by  the  Newton-Raphson  method  now  requires  the  solution 
of  the  second-truncation  equations  an  additional  seven  times  for  each 
iteration.  Clearly,  if  many  iterations  are  required  to  obtain  the 
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correct  solution,  this  procedure  will  be  quite  time  consuming  and  hence 

expensive.  Thus  It  Is  Important  to  have  an  accurate  guess  for  the  seven 

Initial  values  In  order  to  reduce  the  number  of  Iterations  that  are 

necessary.  The  solution  to  the  f lrst- truncation  problem  provides 

reasonable  values  for  u^,  vq,  pQ,  p2  and  Tq.  In  order  to  obtain  a 

first  approximation  to  u^  (0),  p^(0)  and  T2  (0),  It  is  advantageous 

n  n 

to  solve  equations  (3.21)  using  the  first- truncation  solution  to  provide 

values  for  u^»  vq,  pQ ,  p^,  and  Tq  .  Except  for  several  terms  In  the 
dP4 

equation  for  ,  these  equations  are  linear  In  u^,  p^,  and  T2  . 

If  these  nonlinear  terms  are  omitted  and  If  the  shock  parameters  Aq, 

A,,  0 i  are  kept  at  the  first-truncation  values,  one  Iteration  on 

equations  (3.21)  converges  due  to  the  linearity  of  the  problem.  This 

preliminary  solution  provides  an  approximation  for  u^  (0),  p^(0)  and 

H 

T2(0)  whose  accuracy  depends  on  the  accuracy  of  the  first  truncation. 

If  the  first-truncation  solution  is  accurate,  this  preliminary  step 

described  above  can  result  in  a  considerable  reduction  In  computation 

time.  It  will  be  seen  in  the  next  section  that  the  accuracy  of  the 

first  truncation  solution  for  a  common  class  of  problems  can  be  greatly 

improved  by  using  a  more  realistic  approximation  for  ^  than  0^  ■  0. 

The  time  required  to  solve  these  equations  numerically  depends  very 

strongly  on  the  Individual  problem.  Problems  for  large  values  of  the 

Reynolds  number  require  more  time  than  those  for  moderate  values  of  Re. 

The  accuracy  of  the  Initial  guess  for  the  unknown  boundary  conditions 

is  also  quite  Important.  Typical  computation  times  for  the  first-trun- 

2 

cation  solution  range  from  15  seconds  per  iteration  at  Re  -  10  to 

D 


50  seconds  per  iteration  at  Re  *10  .  For  the  second -truncation  problem 

S 

2 

these  times  are  30  seconds  per  iteration  at  Re  -  10  and  200  seconds 

s 

4 

per  Iteration  at  Re  ■  10  .  The  number  of  iterations  that  are  necessary  is 

B 

usually  small  (two  or  three)  unless  the  guessed  values  for  the  initial 
data  are  very  poor.  The  large  computation  times  at  large  values  of 
Re  are  a  consequence  of  having  to  maintain  a  small  step  size  in  the 

S 

integration  procedure  while  integrating  across  the  "inviscid"  region  of 

the  shock  layer,  despite  the  fact  that  the  unknowns  do  not  vary  rapidly 

in  that  region.  This  feature  appears  to  be  a  consequence  of  solving 

for  the  second  derivatives  u  and  T  as  principle  terms.  In  the 

nn  nn  r 

outer,  "lnvlscld1*  region,  u  and  T  must  be  considerably  smaller  than 

the  other  terms  (in  boundary -layer  theory,  they  go  to  zero  exponentially 

as  the  edge  of  the  boundary  layer  is  approached).  If  a  large  number  of 

computations  are  to  be  made  at  large  values  of  Re  ,  it  may  he  worth- 

s 

while  to  devise  a  scheme  that  will  reduce  the  equations  (3.4)  and  (3.5) 
to  first  order  in  the  outer  flow. 

A  possible  reduction  in  computation  time  may  be  achieved  by  origi¬ 
nating  the  integration  of  the  differential  equations  at  the  shock  instead 
of  at  the  body  surface.  For  the  first-truncation  problem  it  would  then 

be  necessary  to  guess  only  three  initial  conditions:  is  ,  u^  (A  )  and 

n 

T  (A  )  .  Thus  an  iteration  based  on  this  alternate  procedure  requires 
n 

one  less  solution  of  the  differential  equations  than  does  the  procedure 
described  earlier.  A  similar  reduction  would  occur  in  the  second-trunca¬ 
tion  problem  since  it  v  ji-i  be  necessary  to  guess  only  six  values  at  the 
shock  in  order  to  initiate  the  integration. 
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D-  Results 


Several  examples  are  considered  in  this  section  in  order  to  evaluate 
the  accuracy  of  the  method  of  solution  and  the  accuracy  of  the  basic  flow 
model.  In  addition,  a  comparison  of  the  first-  and  second-truncation 
results  provides  a  test  of  the  validity  of  local  similarity. 

1.  A  Comparison  to  Second-Order  Boundary -Layer  Theory. 


The  flow  around  a  sphere  at  M  =10,  Y  =  1*4,  a  =  0-7  and 

00 

l/2 

H  =  T  is  used  to  test  the  accuracy  of  the  series -truncation  method- 
For  this  example,  the  body  temperature  is  constant  at  0.6  of  the 
inviscid  stagnation  temperature;  i.e.. 


b  =0.6  and  b„  =  0 
0  2 


where  b  = 


=  b  +  b_sin  s  +■ 
o  2 


(3.22) 


The  results  of  a  boundary-layer  analysis  by  Davis  and  Fliigge-Lotz 
[12]  are  usee*  as  a  standard  of  comparison  for  the  series -truncation 
results.  These  boundary -layer  computations  were  made  for  both  the 
first-  and  second -order  boundary- Layer  equations  by  use  of  an  implicit 
finite-difference  scheme.  A  numerical  scheme  of  this  nature  has  been 
shown  to  provide  accurate  solutions  to  the  boundary-layer  equations . 

The  inviscid  pressure  distribution  on  which  the  boundary-layer  computa¬ 
tions  were  based  was  provided  by  H.  Lomax  of  the  Ames  Research  Center 
of  NASA.  The  method  of  solution  for  the  inviscid  flow  is  described  in 
reference  20  and  has  been  shown  to  be  quite  accurate  [25].  Thus  the 
solution  of  the  first-  and  second-order  boundary -layer  equations 
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computed  by  Davis  and  Flugge-Lots  may  be  considered  to  be  exact  and  to 
provide  a  standard  against  which  the  series-truncation  method  may  be 
measured. 

The  results  for  the  shear  stress,  x/e  {see  equations  (2.22) 

and  (3»l4)),  are  compared  in  figure  3.1.  In  order  to  have  a  significant 

difference  between  the  first-  and  second-order  boundary-layer  results, 

the  rather  low  value  of  100  was  chosen  for  the  shock  Reynolds  number, 
p*U*a* 

Re  =  —  .  Note  that  the  second-order  boundary-layer  effects 

^sh 

decrease  the  wall  shear  slightly  for  this  example.  Since  the  truncation 
results  are  in  the  form  of  a  power  series  centered  at  the  axis,  s  =  0, 
the  results  can  be  expected  to  agree  only  near  s  =  0.  However,  from 
figure  3*1  it  can  be  seen  that  the  first  truncation  gives  a  poor  result 
even  at  the  axis:  the  curve  that  represents  the  first -truncation  result 
has  an  incorrect  slope  at  s  =  0.  The  second-truncation  result  appear* 
to  have  the  correct  slope  at  s  =  0,  but  the  agreement  does  not  extend 
over  any  significant  distance.  The  second  truncation  does,  however, 
give  a  far  more  accurate  solution  than  does  the  first  truncation. 

The  shear-stress  results  can  be  more  precisely  evaluated  if  they 
are  compared  to  a  series  solution  of  the  boundary- layer  equations  rather 
than  to  the  finite -difference  solution  shown  in  figure  3-1.  The 
expression 

x{s,0)/ €  =  1.2  sin  s  -  0.6  sin^s  (3-23) 

agrees  to  0{s^)  with  the  two-term  series  solution  of  the  first-order 
boundary- layer  equations  presented  by  Davis  and  Fliigge-Lotz.  The 
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‘2™  ORDER  BOUNDARY  LAYER 


$,  DISTANCE  FROM  AXIS 

Figure  3.1 

Shear  stress  on  a  sphere  at  Y  *  1.4,  M  =  10,  Re"*  =  100, 

oo  S 

b  =  0.6,  a  =  0.7,  and.  co  =  l/ 2  (e  =  0.1l8);  comparison  of 
truncation  results  with  the  boundary-layer  results  of  Davis 
and  Flugge-Lotz  [12]. 


[•>] 


second-order  boundary-layer  effects  would  decrease  the  coefficients  of 
this  expression  slightly.  However/  in  figure  3.1  these  effects  appear  to 
be  smaller  than  the  computational  errors  associated  with  the  truncation 
method  and  are  ignored  for  the  moment.  The  present  investigation  yields 

T^fye  =  1.50  sin  s  (3- 24a) 

and 

T^tye  =  1.29  sin  s  -  0.51  sin^s  (3 -24b) 

for  the  first  and  second  truncations  respectively.  Thus  despite  the 
pronounced  improvement  of  the  results  shown  by  the  second  truncation/ 
there  is  still  an  inaccuracy  in  the  value  of  t^.  However/  it  should 
be  noted  that  these  errors  are  somewhat  less  than  is  indicated  in 
figure  3*1:  the  curve  representing  the  series  solution  to  the  boundary- 
layer  equations  (equation  3*23)  lies  above  the  finite -difference  curve 
shown  in  figure  3.1  just  as  the  series  truncations  do. 

In  contrast  to  the  wall-shear  results/  the  results  shown  in 
figure  3 *2/  for  the  heat-transfer  rate,  q/e  (see  equations  (2.23) 
and  (3*15))/  agree  with  the  boundary-layer  values  remarkably  well. 

The  first-truncation  value  is  about  five  percent  too  high  at  s  =  0. 


^  The  fact  that  the  f irst-truncation  result  for  q  is  a  constant  is  a 

consequence  of  having  expanded  the  temperature  into  T(s,r?)  = 

2 

T  (n)  +  T9(n)sin  s  +...  .  If  an  expansion  such  as  T(s/n)  = 

T  (njeos^s  +  T„(n)sin  s+.«.  had  been  used,  the  first-truncation 

®  ^  (l)  €  ( ’ )  2 

problem  would  remain  the  same  but  would  yield  ql  '  =  —  q^'ccs  s. 

This  function  would  then  agree  reasonably  well  with  the  boundary- 

layer  values  over  a  more  significant  range  of  values  of  s . 
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This  small  inaccuracy  in  qQ  is  corrected  by  the  second-truncation 

computation.  In  addition*  the  second-truncation  heat-transfer  rate* 

(2)  (2)  2 

‘Jq  +  'sin  s*  agrees  extremely  well  with  the  second-order  boundary- 
layer  result  over  the  entire  s -interval  for  which  the  boundary-layer 


result  is  available. 


The  relative  inaccuracy  of  the  shear-stress  values  seems  to  be 

strongly  related  to  the  description  of  the  shock  6hape.  An  examination 

of  equations  (3*8)  shows  that  a  change  in  the  shock  angle*  0  ^*  affects 

the  value  of  u^  at  the  shock*  i*e.,  u^A^)*  but  not  the  values  of 

T  (A  )  and  p  (A  ).  The  first -truncation  solution  of  this  problem  is 
o  0  o  o 

based  on  an  assumed  value  of  zero  for  0^,  and  it  leads  to  a  computed 

value  of  I.085  for  the  shock  standoff  distance  A  .  The  se  jnd- 

o 

truncation  solution  yields  calculated  values  for  both  0  ^  and  Aq: 


A  =  1.232  and  0  .  =  0.1236 
o  1 


The  effect  of  these  changes  in  the  shock  shape  can  be  evaluated  from 
figure  3*3  where  the  first-  and  second-truncation  values  of  the  flow 
variables  u^*  p^*  and  are  shown.  It  can  be  seen  that  the  changes 

in  0  ^  and  A^  influence  the  value  of  u^  across  the  entire  shock 
layer,  but  near  the  body  surface  the  temperature  and  pressure  profiles 
are  affected  only  slightly  by  the  change  in  value  of  Aq.  For  a  very 
cold  wall,  it  may  be  expected  to  have  a  slightly  larger  effect  on  the 
temperature  than  is  shown  in  figure  3»3«  However,  from  this  example, 
it  is  evident  that  the  ma^or  errors  of  the  first  truncation  should  be 
expected  to  occur  in  the  tangential-velocity  profile  and  hence  in  the 


DISTANCE  FROM  BODY 


Figure  3*3 

Flow  variables  near  the  stagnation  streamline  of  a 
sphere  at  y  =  1.4,  =  10,  Reg  =  100,  b  =  0.6, 

a  =  0.7,  and  oo  =  l/2. 
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wall  shear.  Further,  these  errors  appear  to  be  primarily  a  result  of 
truncating  the  series  that  describes  the  shock  angle  0 . 

This  last  assertion  requires  further  verification,  however,  since 
the  truncation  of  the  "excess"  variables  from  the  normal  momentum  equa¬ 
tion  is  also  a  source  of  error.  As  a  check,  the  first -truncation  equations 
have  been  re-solved  for  this  example.  However,  instead  of  setting 
0^=0,  the  value  found  from  the  second-truncation  solution  is  used; 
i.e.,  a  value  of  0^  =  0.1236  is  specified.  The  resulting  tangential- 
velocity  profile  is  shown  in  figure  3-3  for  comparison.  This  profile 
corresponds  to  a  value  of  1.22  for  The  temperature  profile  and 

the  heat-transfer  rate  do  not  differ  appreciably  from  the  values  pre¬ 
viously  given.  It  can  be  concluded  that  the  major  error  in  the  first 
truncation  originates  from  the  approximation  for  0  the  influence  of 
the  truncation  of  and  v^  from  the  normal  momentum  equation, 

(3*5c),  though  not  negligible,  is  considerably  smaller. 

The  influence  of  the  shock  shape  on  the  second-truncation  problem 

is  somewhat  more  complex.  As  noted  in  section  C,  the  second-truncation 

problem  can  be  made  determinate  in  several  ways.  The  second -truncation 

results  that  were  presented  in  figures  3.1-3-3  and  in  equation  (3.24b) 

were  based  on  an  assumption  that  0^  =  =  0.  In  figure  3-4  the 

effect  of  assuming  only  the  value  of  A^  and  computing  the  consistent 

value  of  0^  from  equation  (3.13)  is  shown.  Values  of  0  , 

(2) 

and  are  plotted  as  functions  of  A^.  For  comparison,  the  value 

(2)  (2) 

of  and  '  given  in  equation  (3.24b)  for  =  0  =  0  are 

also  shown.  As  expected,  the  value  of  t  is  strongly  influenced  by 
the  choice  of  A^  since  u^(A^)  is  a  function  of  0^  and  thus  a 


function  of  A^  (see  equation  (J.lla)).  More  surprising  is  the  noticable 
influence  of  A^  on  T^.  However,  values  of  A^  greater  than  0.04 
begin  to  have  a  noticable  effect  on  the  heat-transfer  results .  Hence 
if  the  outstanding  heat-transfer  results  were  not  simply  fortuitous,  the 
shock  errors  alone  do  not  account  for  all  the  difference  between  the 
boundary -layer  results  and  the  second-truncation  solution. 

It  should  be  noted  that  at  low  Reynolds  numbers  the  present  formula¬ 
tion  of  the  problem  cannot  be  expected  to  yield  results  that  are  identical 
to  those  of  the  second-order  boundary-layer  theory.  Although  the  equa¬ 
tions  of  Chapter  II  are  uniformly  valid  to  0(e)  throughout  the  shock 
layer  and  thur  contain  all  second-order  boundary-layer  effects,  the 
current  formulation  does  not  isolate  these  effects  and  analyze  each  one 
separately  as  is  done  in  boundary-layer  theory.  Consequently,  the 
present  formulation  contains  some  higher-order  effects  that  are  not 

contained  in  the  second-order  boundary-layer  results.  These  differences 

2 

in  the  two  formulations  of  the  problem  are  0(e  )  and  can  be  expected 

to  diminish  as  the  Reynolds  number  is  increased. 

The  variation  with  Reynolds  number  of  the  ratio  of  the  wall  shear 

to  the  first-order  boundary-layer  wall  shear  is  shown  in  figure  3-5  • 

It  can  be  seen  that  the  second-truncation  results  do  indeed  approach  the 

boundary- layer  results  as  the  Reynolds  number  is  increased,  but  somewhat 

slowly.  At  a  Reynolds  number  of  10^,  e  =  0.037  and  the  difference 

-3  2 

between  the  two  formulations  should  be  0(10  ).  At  Res  =  10  ,  the 

.2 

difference  should  be  0(10  ).  Yet,  in  both  cases,  the  actual  differences 
in  the  shear-stress  values  are  about  an  order  of  magnitude  larger  than 


Figure  3.5 

Variation  of  the  stagnation-point  shear  stress  with  the  shock  Reynolds  number 


expected.  Note,  however,  that  the  difference  is  only  a  few  percent  at 

the  larger  values  of  Re  and  is  within  the  variation  due  to  the  un- 

s 

certainty  about  the  shock  shape. 

Accepting  the  accuracy  of  the  results  of  Davis  and  Flugge-Lotz,  we 
arrive  at  the  following  conclusions. 

a.  The  first  truncation  results  are  considerably  less  accurate 
than  those  of  the  second  truncation,  even  at  the  axis  of 
symmetry. 

b.  The  heat -transfer  rate  is  determined  quite  accurately,  but 
significant  errors  occur  in  the  values  of  shear  stress. 

c.  The  error  in  the  shear-stress  values  is  primarily  due  to  the 
sensitivity  of  the  tangential  velocity  to  approximations  for 
the  shock  shape. 

d.  The  errors  due  to  the  truncations  increase  as  the  Reynolds 
number  is  decreased. 

2.  Local  Similarity 

The  results  of  the  first  example,  particularly  those  shown  in 
figure  3-5,  are  of  interest  for  an  evaluation  of  the  results  of  local- 
similarity  analyses,  e.g.  Ho  and  Probstein  [19],  Probstein  snd  Kemp 
[28],  and  Goldberg  [17] •  In  each  of  these  investigations,  comparisons 
to  boundary -layer  results  much  like  those  in  figure  3-5  are  presented. 
In  each  case  the  low-Reynolds-number  results  are  qualitatively  very 
similar  to  those  presented  above;  that  is,  the  shear-stress  values  are 
found  to  be  considerably  larger  than  predicted  by  first-order  boundary- 
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l&yer  theory.  Further,  the  method  of  analysis  that  is  used  in  these 
investigations  is  very  similar  to  the  first -truncation  analysis.  The 
variables  have  been  assumed  to  be  "locally  similar"  in  the  region  of 
the  axis  and  hence  are  assumed  to  be  adequately  described  by  expressions 
such  as  u  =  s  u-^n),  v  =  {l  -  2  s  )  vQ(n),  p  =  {l-s  )  pQ(n)  -  ^  s  p2(n), 
etc.  The  similarity  to  the  expressions  used  in  the  first  truncation  is 
evident.  In  addition,  these  investigations  have  assumed  that  the  body 
and  shock  are  concentric  {0^  =  0),  a  fact  that  has  been  shown  to  be  of 
considerable  importance.  The  results  of  the  first  example  of  this 
section  imply  that  a  portion  of  the  divergence  from  boundary -layer  theory 
found  in  these  investigations  is  a  consequence  of  the  method  of  analysis 
and  not  a  property  of  the  low -Reynolds -number  flow.  That  is,  the  use  of 
the  method  of  local  similarity  leads  to  errors  which  exaggerate  the 
differences  with  the  results  of  boundary-layer  theory. 

A  verification  of  the  above  conclusions  by  a  direct  comparison  of 
the  results  of  a  series -truncation  analysis  to  those  of  the  local- 
similarity  analyses  is  somewhat  difficult.  The  series-truncation 
analysis  has  been  applied  to  the  problem  defined  by  Ho  and  Probstein, 
i-e.,  Y  =  — t  =  °°,  a  ~  0-71*  b0  =  0.05,  b2  =  0.0,  =  T^^,  and 

no-slip  boundary  conditions  at  the  body  surface.  However,  Goldberg 
gives  only  a  range  of  values  cf  flight  velocities,  altitudes,  and  body 
temperatures  at  which  his  computations  were  made.  Thus  there  is  some 
uncertainty  about  what  values  of  Mach  number,  Prandtl  number,  and  body 
temperature  should  oe  used  in  the  series -truncation  analysis  for 
comparison  with  his  results.  However,  Goldberg  does  present  a  result 
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for  - —  =  0.10,  which  corresponds  to  the  choice  of  Y  =  -jr  above, 

Psh  y 

and  he  specifies  no-slip  boundary  conditions.  Further,  for  large  values 

of  the  Mach  number,  the  precise  value  of  should  have  only  a  minor 

effect;  and  the  normalization  of  the  results  with  respect  to  first- 

order  boundary -layer  values,  as  wan  done  by  Goldberg,  should  minimize 

the  influence  of  the  unknown  parameters.  Hence  a  comparison  of  the 

shear-stress  results  of  Goldberg,  Ho  and  Probstein,  and  the  series - 

truncation  method  is  given  in  figure  3.6.  The  results  of  Ho  and 

Probstein  and  of  the  series  truncation  have  been  normalized  by  t  = 

l.liyReg  sin  s  +•••  .  This  function  was  obtained  from  reference  19 

and  is  based  on  an  inviscid  pressure  gradient  that  was  obtained  from  a 

thin-shock-layer  analysis.  Interpretation  of  these  results  is  complicated 

by  the  general  lack  of  agreement  among  the  several  results. 

Goldberg's  results  illustrate  two  effects  not  contained  in  the 

P» 

other  results.  First,  as  the  shock  density  ratio,  - —  ,  becomes 

Psh 

smaller  (requiring  y  ->l),  there  is  a  greater  deviation  from  the 

Peo  5 

first-order  boundary -layer  values  (note  that  the  choice  of  (■—— )  Re^ 

sh 

as  the  abscissa  in  figure  3-6  reduces  the  appearance  of  this  effect). 

This  is  to  be  expected  since  the  shock  layer  becomes  thinner  and  the 
second-order  displacement -thickness  effect  becomes  more  pronounced  than 
it  was  in  the  first  example  of  this  section.  Second,  Goldberg's  results 
show  a  decrease  in  the  shear  stress  at  very  low  Reynolds  numbers.  This 
decrease  is  a  consequence  of  having  used  a  modified  set  of  shock  condi¬ 
tions  which  account  for  the  thickening  of  the  shock.  The  truncation 
results  and  those  of  Ho  and  Probstein  are  based  on  the  Rankine-Hugoniot 
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shock  relations  and  thus  do  not  exhibit  this  effect. 

The  agreement  between  the  first  truncation  and  Goldberg's  result 
P» 

for  =0.10  is  actually  fairly  good  considering  the  uncertainty 

Psh  k 

about  some  of  the  flow  conditions.  At  Re  =  10  ,  Goldberg's  results 

8 

(when  slightly  extapolated)  predict  a  shear  stress  which  is  about  25 

percent  higher  than  the  boundary-layer  value*  and  the  first -truncation 

solution  is  about  17  percent  higher.  In  addition,  computations  have 

shown  that  approximating  factors  of  1+en  by  unity,  as  done  by  Goldberg, 

reduces  the  difference  between  these  two  results  by  about  one  fourth. 

The  second-truncation  results  are  qualitatively  the  same  as  they 

were  in  the  first  example.  The  decrease  in  wall  shear  below  the  first- 

truncation  values  is  somewhat  less  than  it  was  before,  but  is  still 

k 

substantial.  In  fact,  at  Re  =  10  ,  the  second-truncation  result  is 

s 

only  six  percent  higher  than  the  boundary -layer  value  (compared  to 
25  and  17  percent  cited  above).  Since  the  displacement-thickness  effect 
is  larger  for  this  example  than  it  was  for  the  first  example  of  this 
section,  the  total  second-order  boundary -layer  effect  can  be  expected 
to  be  positive,  unlike  the  result  shown  in  figure  3.5.  As  anticipated, 
the  second-truncation  results  provide  a  much  better  agreement  with  the 
boundary- layer  values  than  the  results  of  a  local-similarity  analysis  do. 

A  comparison  to  the  results  of  He  and  Probstein,  however,  is  not 
conclusive  since  the  results  of  the  first  truncation  and  -hose  of  Ho 
and  Probstein  are  not  in  agreement.  The  reasons  for  this  difference 
are  not  fully  understood.  Approximately  one  third  of  the  difference 


(at  Re  =  100)  is  the  result  of  the  omission  of  the  viscous  terms 
8 
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cku'u  T  and  ep(/<  +  j  - )u  from  the  tangential  momentum  equation 

n  r  n 

(see  equation  (2.3))  by  Ho  and  Prcbstein.  On  the  other  hand;  they 
retain  several  viscous  terms  in  the  normal  momentum  and  energy  equations 
which  are  not  included  in  equations  (2.4)  or  (2.5).  The  influence  of  all 
these  terms  except  those  that  are  proportional  to  (^vn)n  ^ias  been 
investigated  and  has  been  found  to  change  the  results  by  less  than  one 
percent.  The  effect  of  (nvn)n  was  not  checked  since  the  computational 
difficulties  described  earlier  in  this  chapter  would  require  extensive 
modification  of  the  computational  procedures.  Thus  the  major  part  of 
the  difference  is  unexplained. 

However,  it  should  be  noted  that  the  pressure  distribution  across 
the  shock  layer  computed  by  Ho  and  Probstein  differs  radically  from  the 
truncation  result  for  pressure.  In  all  the  examples  hat  have  been 
computed,  the  inviscid  pressure  mechanisms  which  are  included  in  equation 
(2.4)  yield  pressure  distributions  like  that  shown  in  figure  (3-3);  i.e., 
the  pressure  increases  monotonically  from  the  shock  to  the  body  surface. 
In  contrast,  the  computations  of  Ho  and  Probstein  yield  a  pressure 
distribution  that  has  a  rather  steep  drop  immediately  behind  the  shock 
before  increasing  to  the  surface  value  (for  an  example,  see  figure  3-7)> 
As  pointed  out  by  Shih  and  Krupp  [30],  this  decrease  in  pressure,  which 
is  apparently  a  result  of  viscous  effects,  is  suspect  because  one  does 
not  expect  the  viscous  forces  to  be  larger  immediately  behind  the  shock 
than  near  the  body  surface.  The  investigation  by  Kao  [21]  has  yielded 

The  term  j  (k  +  j  sin  S/r)Tn  is  also  omitted  from  the  energy 
equation.  This  term  has  a  significant  influence  only  on  the  heat 
transfer,  however* 
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a  similar  pressure  distribution,  but  Levinsky  and  Yoshihara  [24]  and 
Goldberg  have  obtained  pressure  distributions  that  do  not  show  this 
effect  despite  their  use  of  equations  that  are  quite  similar  to  those 
used  by  Ho  and  Probstein. 

It  may  be  possible,  of  course,  that  the  method  of  local  similarity 
yields  a  more  accurate  solution  to  the  equations  of  reference  19  than 
it  does  to  those  of  Chapter  II.  However,  it  seems  more  likely  that  a 
series -truncation  analysis  <  “  the  equations  used  by  Ho  and  Probstein 
would  yield  results  like  those  exhibited  in  the  two  examples  considered 
in  this  section;  i.e.,  the  local  similarity,  or  first  truncation, 
results  would  be  significantly  improved  a.t  the  axis  by  the  more  complete 
analysis  of  the  second  truncation.  This  conclusion  is  further  sub¬ 
stantiated  by  the  following  discussion  of  the  analysis  by  Kao:  Kao's 
results  are  similar  to  those  of  Ho  and  Probstein,  and  it  will  be  shown 
that,  in  general,  his  analysis  should  yield  substantial  differences 
between  the  first-  and  sec end -truncation  results. 

3 •  Influence  of  the  Body  Temperature 

The  results  that  have  been  described  above  appear  to  be  in  contradic¬ 
tion  to  the  results  of  an  investigation  by  H.C.  Kao  [21] .  He  solved 
the  Navier-Stokes  equations  by  use  of  the  method  of  series  truncation 
and  fiund  very  little  difference  between  the  results  of  the  first  and 
second  truncations.  Hence  he  concluded  that  the  method  of  local  similarity 
was  valid  (in  the  sense  that  it  y  ,us  correct  values  for  the  first 
term  of  a  more  complete  series  so]  „ion).  The  contradiction,  however, 
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is  more  apparent  than  real  and  can  be  resolved  by  a  comparison  of  the 
computational  procedure  of  Kao  to  that  used  in  this  chapter. 

The  truncation  procedure  used  in  this  chapter  has  been  adapted 
from  Kao's  analysis,  and  hence  the  two  formulations  of  the  problem 
are  quite  similar.  However,  since  Kao  solved  the  inverse  blunt-body 
problem,  the  location  of  the  body  surface  appears  as  an  unknown  quantity 
instead  of  the  shock  location.  In  the  integration  of  the  equations 
inward  from  the  shock,  the  location  of  the  body  can  be  determined  from 
one  of  the  wall  boundary  conditions.  However,  Kao  adopted  an  alterna¬ 
tive  procedure  of  specifying  the  body  location  in  which  case  the  surface 
temperature  cannot  be  imposed  as  a  boundary  condition  but  must  be 
considered  to  be  a  result,  much  like  the  shear  stress  and  heat  transfer. 

More  specifically,  if  both  T,  (see  equation  (3. 7))  and  the  body 

b2 

location  away  from  the  axis  of  symmetry  are  specified  in  the  second- 
truncation  problem,  the  equations  are  over-determined.  Hence,  Kao 
specifies  that  the  body  and  shock  be  concentric  and  obtains,  as  part 
of  his  solution,  the  temperature  distribution  of  the  body  surface.  Thus 
Kao's  conclusion  about  local  similarity  is  valid  only  if  the  tempera¬ 
ture  of  the  body  is  such  that  the  body  and  the  shock  are  concentric. 

The  temperature  distribution  which  produces  this  spherical  symmetry 
is  one  in  which  the  body  is  strongly  cooled  downstream  of  the  axis  of 
symmetry.  Kao  solved  the  problem  defined  by  Y  =  >  M  =  10, 

y  « 

Re  =  10,  <J  =  0.7>  b  -  0.048,  and  no-slip  boundary  conditions.  The 
resulting  temperature  distribution  is  given  by  bg  =  -0.5l6.  When  the 
present  series -truncation  analysis  is  applied  to  this  problem  (with  bg 
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Kao 

FIRST  and  SECOND  TRUNCATION 


Figure  3.7 

Comparison  with  the  flow  variables  computed  by  Kao  [21]; 

Y  =  ll/9,  M  =  10,  Re  =  10,  b  =  0.048,  b„  =  -0.516, 

00  8  O  c 

o' =  0.7,  o)=l/2>  and  no-slip  boundary  conditions  (e  =  0.369)* 


specified  as  -0.516),  It  is  found  that  the  equations  of  Chapter  II 

yield  a  nearly  spherically  symmetric  problem  for  this  temperature  dis- 

A(2> 

tribution  (more  specifically,  we  obtain  =  1.003  and  0  ^  =  O.OO36). 

Ao 

Despite  this  agreement  with  the  results  of  Kao,  the  flow  variables  ob¬ 
tained  from  these  two  solutions  do  not  agree  well,  as  seen  in  figure 
3 *7*  The  choice  of  Reg  =  10  is  undoubtedly  outside  the  range  of 
validity  of  the  basic  flow  model  of  this  investigation  (note  that  the 
viscous  effects  extend  across  the  entire  shock  layer  ).  Further,  the 
radically  different  pressure  distributions  that  were  discussed  with 
regard  to  the  previous  example  are  also  encountered  here  and  preclude 
the  possibility  of  making  more  than  a  brief  comparison. 

The  influence  of  the  downstream  temperature  distribution  has  also 
been  Investigated  for  the  flow  problem  defined  in  the  first  example; 

i.e.,  Y  =  1.4,  M  =  10,  Re  =  100,  cr  =  0.7,  and  b  =  0.6.  In 
00  s  o 

figure  3*8a,  the  variation  of  and  q^  with  is  shown.  As  is 
to  be  expected,  the  values  of  and  are  strongly  dependent  upon 

the  value  of  b^.  In  addition,  however,  the  results  of  these  computa¬ 
tions  clearly  show  the  presence  of  an  upstream  influence  in  the  flow 
since  t^,  q^,  and  Aq  are  significantly  influenced  by  the  downstream 
wall  temperature.  These  quantities  are  shown  in  figure  3.8b  where  they 
have  been  normalized  by  their  respective  values  computed  in  the  first 
truncation.  Also  shown  is  the  value  of  0^.  In  the  first -truncation 
computation,  the  value  of  b^  is  of  no  concern.  If  no  coupling  occurred 
between  different -ordered  coefficients  of  the  series  expansions 
(equations  (3*1))/  the  first-truncation  results  would  be  valid  for  all 
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Figure  3.8a 

Influence  of  the  downstream  body  temperature  (at  Y  =  1.4, 
Mm  =  10,  Rea  =  100,  cr  =  0.7>  and  tu  =  l/2)on  the  second 
coefficients  of  the  expansions  for  the  shear  stress  and 
heat -transfer  rate. 
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0  -IjO  -2.0  -3.0  -4.0  -5.0 


Figure  3.8b 

Influence  of  the  downstream  body  temperature  (at  Y  =  1.4, 

Mw  =  10,  Reg  =  100,  u  =  c.7,  and  <0  =  l/2)  on  the  shock  shape, 
and  the  first  coefficients  of  the  expansions  for  the  shear 
stress  and  heat -transfer  rate. 
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values  of  Hence  the  deviation  of  these  normalized  quantities  from 

a  value  of  1  represents  an  error  in  the  first  truncation  due  to  an 
upstream  influence  in  the  flow.  It  can  be  seen  that  the  magnitude  of 
thi6  error  depends  on  the  value  of  the  downstream  temperature,  i.e., 
on  b,>.  For  the  case  of  b^  ~  -3-8,  the  body  and  shock  are  essentially 
concentric  since  0  ~  0 .  For  this  case,  our  computation  is  equivalent 

to  the  computational  method  used  by  Kao.  It  can  be  seen  that  the  changes 
in  shear  stress,  heat  transfer,  and  shock  standoff  distance  are  con¬ 
siderably  smaller  than  they  are  for  the  constant  wall-temperature  case. 
The  smallest  changes  occur  at  slightly  larger  values  for  b  and  0  , 
but  it  is  obvious  that  the  overall  agreement  between  the  first-  and 
second-truncation  results  is  considerably  improved  by  requiring  the 
shock  and  body  to  be  concentric.  Hence,  as  before,  we  find  that  there 
is  a  temperature  distribution  that  yields  a  spherically  symmetric  flow, 
and  for  this  case,  the  local-similarity  analysis  is  fairly  accurate. 
However,  the  required  temperature  distribution  differs  quite  markedly 
from  the  constant  wall  temperature  which  is  commonly  used  as  a  boundary 
condition. 

It  has  also  been  found  that- the  specification  of  an  adiabatic  wall 

condition  leads  to  results  much  like  those  obtained  from  a  constant  wall 

temperature.  Specifying  adiacatic  conditions  with  the  flow  parameters 

of  the  above  example  leads  to  a  temperature  distribution  of  b  =  O.989  - 
o 

0.100  sin  s.  The  magnitude  of  the  changes  between  the  first  and  second 

A<2>  r<2> 


truncations  are  described  by  =  1-138,  0^  =  0.128,  and  "yyy  =  °-851. 


In  simmary,  we  conclude  that  the  results  of  the  method  of  local 
similarity  are  accurate  if  the  geometry  of  the  flow  boundaries  has  a 
spherical  symmetry.  However,  the  most  commonly  used  boundary  conditions 
on  the  temperature,  an  adiabatic  wall  and  a  constant  wall  temperature, 
do  not  correspond  to  flows  having  such  symmetry. 

4.  Influence  of  the  Shock  Thickness 

The  simplification  of  the  Navier-Stokes  equations  made  in  Chapter 
II  and  the  use  of  the  Rankine-Hugoniot  Bhock  relations  restrict  the 
accuracy  of  the  flow  model  at  very  low  Reynolds  numbers.  A  comparison 
to  the  results  of  Levinsky  and  Yoshihara  [24]  provides  some  insight 
into  the  range  of  validity  of  the  flow  model  used  in  this  investigation. 
The  equations  that  are  used  in  reference  24  are  also  a  simplified  form 
of  the  Navier-Stokes  equations,  similar  to  those  used  in  references  19 
and  28.  Since  these  equations  omit  several  terms  of  0(c)  (see  page  So), 
they  cannot  be  expected  to  test  the  validity  of  the  omission  of  terms  of 
0(c  )  from  equations  (2.2)-(2.7)«  However,  more  important  is  the  fact 
that  no  discontinuous  shock  is  used  in  the  formulation  of  the  flow 
model  by  Levinsky  and  Yoshihara.  Instead,  the  differential  equations 
are  applied  from  the  body  to  the  free  stream,  and  the  effect  of  a 
thickening  shock  at  low  Reynolds  numbers  is  obtained.  Although  we 
have  noted  that  the  method  of  solution  which  is  used  —  the  method  of 
local  similarity  --  may  not  be  particularly  accurate,  we  have  seen  that 
it  is  equivalent  to  the  first-truncation  analysis,  and  therefore 
meaningful  comparisons  can  be  made.  In  addition,  an  inability  to 


adequately  describe  the  shock  shape  has  been  shown  to  he  responsible  for 

much  of  the  inaccuracy  of  the  first  truncation.  Hence,  it  is  of  interest 

to  investigate  whether  similar  inaccuracies  occur  when  the  flow  model 

docs  not  contain  a  discontinuous  shock. 

The  problem  considered  by  Levinsky  and  Ycshihara  is  defined  by 
5  l/2 

y  =  7,  M  =10,  ct=0-75j  p=T  ,  and  no -slip  boundary  conditions. 

j  <*> 

Both  cold  walls  and  adiabatic  walls  were  considered.  Three  values  of  the 
free-stream  Reynolds  numbers  were  used:  13,652;  1,302;  and  152.  These 
values  correspond  to  a  shock  Reynolds  number  of  2409,  244,  and  26.8, 
respectively.  The  first -truncation  analysis  of  this  chapter  has  been 
applied  to  the  adiabatic  wall  case,  and  the  results  are  shown  in  figures 
3-9  a-c.  The  functions  shown  in  these  figures  were  defined  by  Levinsky 
and  Yoshihara  and  are  given  in  terms  of  the  truncation  variables  by 


_  p0  -  p0+p2 

ui '  T  -  eyo  ’  P  -  7-  ’  h  -  -ITT 

AW  AW 


where  p  i6  the  stagnation-point  density  for  adiabatic  wall  condi- 

An 

tions.  The  last  relation  shove  is  a  consequence  of  the  pressure  being 
represented  by  an  expression  of  the  form  p(s,n)  = 

PAy[p(n)c°8  s  +  p2(n)sin  sj  in  reference  24  (compare  to  equation 
(3.1c)). 

At  the  largest  value  of  Reynolds  number  considered,  Re  =  2409, 

s 

the  results  of  Levinsky  and  Yoshihara  exhibit,  a  very  thin  shock- 
transition  region  and  a  distinct  separation  of  a  narrow,  viscous  boundary 
layer  from  the  main  region  of  inviscid  flow  (figure  3*9a).  The  first- 
truncation  results  agree  reasonably  well  except  for  the  location  of 
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Figure  J.9a 

Comparison  with  the  results  of  Levinsky  and  Yoshihara  [24] 

at  Y  =  5/3^  =  10,  <r  -  0.75,  tu  =  l/2,  adiabatic  wall, 

no-slip  boundary  conditions,  and  Re  =  2409  (e  =  0.024). 
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the  shock.  This  discrepancy  has  been  found  to  be  a  result  of  an  approx¬ 
imation,  1  +  en  ~  1,  which  was  used  in  the  equations  of  reference  24. 
When  this  approximation  is  made  in  the  first-truncation  problem,  the 
shock  position  corresponds  to  the  outer  edge  of  the  shock-transition 
region  obtained  by  Levinsky  and  Yoshihara.  This  has  been  illustrated 
in  figure  3*9a  by  the  results  for  the  pressure.  Hence,  as  is  to  be 
expected,  the  agreement  is  quite  good  at  large  Reynolds  numbers  since 
the  shock  thickness  is  quite  small.  Unless  otherwise  noted,  the  series - 
truncation  computations  that  follow  do  not  use  the  approximation 
1  +  en  ~  1. 

At  Re  =  244,  the  shock  has  thickened  considerably,  and  the 
8 

"boundary-layer"  is  a  substantial  portion  of  the  shock  layer,  as  shown 
in  figure  3* 9b-  The  variables  still  match  rather  well,  however,  except 
for  the  details  of  the  shock-transition  region;  the  effects  of  the 
thickening  of  the  shock  do  not  appreciably  alter  the  flow  in  the  viscous 
layer.  The  results  cited  earlier  in  this  section  have  shown  that  the 
assumption  of  f*  ^  =  0  in  the  first  truncation  is  quite  influential. 

The  good  agreement  shown  in  figure  3*9b  indicates  that  an  analogous 
assumption  must  be  inherent  to  the  analysis  by  Levinsky  and  YoshVnara 
despite  the  fact  that  there  is  no  discontinuous  shock  in  their  flow 
model.  The  second-truncation  solution  is  again  significantly  different 
from  the  first-truncation  solution,  as  is  shown  in  figure  3*9b  by  the 
result  for  the  tangential  velocity. 

At  Re  -  26.8,  the  shock-transition  region  and  the  "boundary  layer" 
have  merged,  as  shown  in  figure  3 .9c;  and,  of  course,  the  mod^l  defined 
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Figure  3*  9b 

Comparison  with  the  results  of  LevinsJiy  and  Yoshihara  (24] 

at  y  =  5/3/  Mm  =  10/  <J  -  0.75/  03  =  l/2,  adiabatic  wall/ 

no-slip  boundary  conditions/  and  Re  =  244  (e  =  0.077). 
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LEV1NSKY  and  YOSHhARA 
FIRST  TRUNCATION 


I 


in  Chapter  II  cannot  describe  this-  The  agreement  with  the  truncation 

results  now  deteriorates,  and  the  influence  of  the  thickened  "shock" 

causes  a  considerable  reduction  in  the  shear  stress  at  the  wall.  (This 

effect  was  previously  noted  in  a  comparison  with  the  results  of  Goldberg 

1 171*)  Hence  at  such  low  values  of  Re  ,  it  is  necessary  to  either 

8 

modify  the  shock  conditions  (e.g.,  references  29,  2,  and  17)  or 

integrate  the  differential  equations  through  the  shock-transition  region. 

In  addition,  it  can  be  seen  that  the  omission  of  terms  of  0(e  )  from 

equations  (2.3)-(2.5)  leads  to  a  not ic able  error  for  such  a  low  Reynolds 

number.  With  no  viscous  terms  in  the  normal  momentum  equation,  (2.5 ), 

dp 

the  pressure  gradient  — ^  is  zero  at  the  wall,  as  in  boundary-layer 
theory.  With  the  addition  of  the  terms 


(nun)g  +p(cot  s)u 


1  +  en 


n( 


to  the  normal  momentum  equation  (as  in  references  17,  19>  21,  24),  the 
normal  pressure  gradient  at  the  wall  becomes 


dn 


n=o 


=  -2e2p(To(0))U;L  (0)  • 


2 


n 

It  can  be  seen  from  figure  3. 9c  that  this  pressure  gradient  becomes 
quite  important  at  this  low  Reynolds  number.  Hence,  it  appears  that 
for  values  of  the  shock  Reynolds  number  of  the  order  of  100  or  lower 


1 

2 


These  terms  are  based  on  an  assumption  of  zerb  bulk  viscosity;. 

This  equation  is  valid  only  for  no-slip  boundary  conditions  at-  the  will. 
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the  formulation  of  the  problem  as  done  in  Chapter  II  can  lead  to  sig¬ 
nificant  errors.  For  these  cases,  the  problem  should  be  formulated  to 

p 

be  valid  to  0(e  ),  including  the  effect  of  a  shock  of  finite  thickness. 

In  passing,  it  should  be  noted  that  the  results  of  Levinsky  and 
Yoshihara  do  not  exhibit  the  type  of  pressure  distribution  across  the 
shock  layer  that  was  obtained  in  references  19  and  21  (and  that  was 
illustrated  in  figure  3*7)  although  the  equations  which  are  solved  are 
nearly  identical  to  those  cf  reference  19* 

E.  Summary 

In  this  chapter,  the  method  of  series  truncation  has  been  used  to 
obtain  solutions  to  the  equations  of  Chapter  II  for  flow  around  shperical 
bodies.  Due  to  the  complexity  of  the  equations,  it  was  not  possible  to 
compute  a  large  number  of  terms  and  to  thus  extend  the  validity  of  the 
solutions  beyond  the  immediate  neighborhood  of  the  axis-  However,  the 
purpose  of  this  investigation  has  been  to  provide  accurate  data  at 
the  axis  to  serve  as  initial  data  for  a  separate  computational  scheme 
(which  is  discussed  in  the  next  chapter).  In  addition,  the  solutions 
have  yielded  considerable  information  about  the  validity  of  the  basic 
flow  model  that  was  adopted  in  Chapter  II  and  about  the  accuracy  of 
the  method  of  local  similarity  —  a  method  that  has  been  used  in 
numerous  investigations  of  the  blunt-body  problem.  The  main  results 
are  summarized  below. 

1.  The  basic  flow  model  that  was  adopted  in  Chapter  II  appears  to 
be  accurate  for  values  of  the  shock  Reynolds  number  down  to  the  order 
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2 

of  100.  Below  that  point,  terms  of  0(e  )  that  have  Been  omitted  from 
the  differential  equations  and  the  effect  of  a  thickening  shock  "become 
increasingly  important. 

2.  The  first -truncation  results  are  not,  in  general,  accurate  even 
near  the  axis  of  symmetry.  The  first  truncation  consistently  yields 

values  of  the  wall  shear  stress  that  are  substantially  too.  large  and 

,  •- 

values  of  the  shock  standoff  distance  that  are  substantially  too  small. 
The  heat-transfer  rate  is  overestimated,  but  the  error  in  this  quantity 
is  moderate. 

The  second-truncation  results  are  considerably  more  accurate  but 
still  contain  noticeable  errors  at  the  axis.  These  errors,  however, 
are  generally  small  for  Reynolds  numbers  within  the  range  of  validity 
of  the  basic  flow  model.  The  main  error  occurs  in  the  shear  stress 
while  the  heat -transfer  results  are  quite  accurate.  The  sensitivity 
of  the  tangential  velocity  component  to  changes  in  the  shock  slope 
and  an  inability  to  adequately  determine  the  shock  shape  are  the  primary 
sources  of  this  error. 

3.  The  large  errors  of  the  first  truncation  indicate  that  the 
method  of  local  similarity  may  also  result  in  quite  large  errors  since 
the  method  of  local  similarity  is  equivalent  to  the  first  approximation 
of  the  series -truncation  analysis. 

4.  The  magnitude  of  the  errors  of  the  first  truncation  are  depen¬ 
dent  on  the  wall-temperature  distribution.  For  a  body  that  is  highly 
cooled  downstream  of  the  axis,  the  body  and  shock  are  concentric,  and 
the  first  truncation  is  very  accurate.  However,  the  most  commonly  used 


boundary  conditions  --  a  constant  wall  temperature  and  an  adiabatic 
wall  --do  not  correspond  to  a  flow  that  has  thiB  spherical  symmetry. 


Chapter  IV 


FINITE  DIFFERENCE  METHODS 

The  investigation  of  the  preceding  chapter  has  provided  solutions 
of  reasonable  accuracy  to  the  equations  of  Chapter  II  but  only  In  the 
vicinity  of  the  axis  of  symmetry.  To  obtain  solutions  that  are  valid 
over  a  larger  region  of  the  shock  layer,  it  is  necessary  to  turn  to 
other  methods.  In  Chapter  II  it  was  noted  that  the  characteristic 
surfaces  of  the  governing  partial  differential  equations  are  real,  and 
hence,  since  the  equations  do  not  have  an  elliptic  character,  they  may 
be  suitable  for  solution  as  an  initial-value  problem.  The  solutions  of 
Chapter  III  would  appear  to  contradict  this  since  they  have  shown  an 
upstream  influence  in  the  flow,  and  such  an  influence  is  generally 
associated  with  an  elliptic  character  in  the  flow  equations.  However, 
for  the  present,  It  is  assumed  that  a  solution  can  be  obtained  by  using 
the  axis  of  symmetry  as  the  initial  line  and  the  solutions  of  Chapter  III 
as  the  initial  data.  A  method  that  Is  ideally  suited  to  the  purpose  of 
finding  such  a  solution  is  one  that  is  based  on  the  use  of  an  implicit 
fin ice -difference  scheme.  Implicit  finite -difference  schemes  have  been 
used  to  provide  accurate  solutions  to  parabolic  partial  differential 
equations  (such  as  the  boundary -layer  equations).  In  addition,  investi¬ 
gations  have  shown  that  the  method  works  quite  well  on  simple  hyperbolic 
equations.  Hence,  it  appears  feasible  to  apply  such  a  method  to  the 
equations  of  Chapter  II. 
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A.  Background 

Two  previous  applications  of  implicit  finite -difference  methods  to 
the  viscous  blunt-body  problem  have  been  made  (other  them  boundary-layer 
analyses).  The  first  of  these  was  by  H.K.  Cheng  (2)  to  obtain  solutions 
to  the  thin-shock- layer  equations.  The  second  application  was  by 
R.  T.  Davis  and  W.  J.  Chyu  (llj  to  the  equations  of  Chapter  II  but  for 
the  case  of  constant  density  throughout  the  shock  layer. 

The  equations  that  were  used  by  Cheng  were  parabolic  in  character, 

and  thus  were  ideally  suited  for  solution  by  the  finite -difference  scheme. 

However,  several  terms  of  0(e)  have  been  omitted  from  the  equations 

of  that  investigation.  Hence,  the  equations  of  Chapter  II  should  provide 

a  more  accurate  representation  of  the  flow  at  moderately  low  Reynolds 

numbers  than  is  provided  by  the  equations  of  reference  2.  The  addition 

of  these  terms  to  the  equations  (and  the  addition  of  slip  and  temperature 

jump  conditions  at  the  wall)  does  not  change  the  parabolic  nature  of  the 

equations  and  thus  should  not  appreciably  add  to  the  difficulty  of 

obtaining  solutions.  However,  two  additional  terms  that  do  not  appear 

in  the  thin-shock-layer  equations  have  been  retained  in  equation  (2.4), 

the  normal  momentum  equation,  and  these  terms  do  have  an  effect  on  the 

puv 

ba^ic  nature  of  the  problem.  These  terms,  pw  and  \rrr  ,  are  of 

n  1+Ken 

order  unity  in  the  inviscid  region  of  the  flow  field  and  thus  are  essential 
to  an  adequate  description  of  the  inviscid  flow  field  (of  course,  this 
conclusion  is  reached  without  regard  to  the  fact  that  the  shock  layer 
is  thin).  A  discussion  of  the  significance  of  these  term6  for  the 
description  of  the  flow  in  the  shock  layer  was  contained  in  Chapter  II. 
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The  addition  of  these  two  terms  to  the  equations  results  in  the  appearance 
of  two  11  non-parabolic"  characteristic  surfaces  that  are  described  by 
equation  (2-9) •  Hence  the  equations  are  no  longer  of  a  purely  parabolic 
type,  and  the  implications  of  this  fact  with  regard  to  the  use  of  a 
finite -difference  scheme  are  unknown. 

In  addition,  the  thin-shock-layer  concept  simplifies  the  description 
of  the  shock  shape  and  the  determination  of  the  shock  position.  The 
location  of  the  shock  can  be  determined  by  a  balance  of  the  mass  flow 
in  the  shock  layer  with  the  mass  flow  that  enters  the  shock  layer  from 
the  free  stream.  Equation  (2.l8)  is  a  formal  representation  of  this 
balance.  The  mass  flow  that  crosses  the  shock  from  the  free  stream 
into  the  shock  layer  is  proportional  to  (r  +  cA  sin  0)^+1.  In  general, 
since  this  expression  is  a  function  of  A,  the  value  of  the  expression 
depends  upon  the  solution  of  the  flow  variables  in  the  shock  layer. 
However,  to  the  first  approximation  of  the  thin-shock-layer  theory, 
the  body  and  shock  coincide,  i.e.  A  -»  0,  and  the  macs  flow  in  the 
shock  layer  is  proportional  to  r^+1.  Thus,  in  the  thin-shock- layer 
approach,  the  locations  of  the  flow  boundaries  are  known  prior  to  solving 
for  the  flow  in  the  shock  layer.  In  addition,  the  shock  angle  is  the 
same  as  the  local  body  angle  (0  =  0),  and  the  values  of  u,  v,  T,  and 
p  at  the  shock  are  also  known  independently  of  the  solution  within  the 
shock  layer.  Thus,  in  the  more  general  case  considered  in  this  chapter, 
the  computation  of  the  boundary  position  and  shape  represents  a  complica¬ 
tion  which  was  not  encountered  by  Cheng:  it  will  be  necessary  to  compute 
the  location  of  the  shock  position  and  slope  at  each  step  since  the 
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amount  of  mass  flow  along  the  shock  layer  ic  not  known  a  priori  as  it 
is  under  the  thin-shock-.layer  assumption.  It  has  already  been  seen  in 
Chapter  II  that  such  a  computation  can  lead  to  difficulties. 


Both  of  these  complications  were  encountered  by  Davis  and  Chyu 
when  they  solved  the  constant -density  case  using  equations  that  were 
equivalent  to  those  in  Chapter  II.  It  was  found  that  valid  results 
could  be  obtained  but  that  two  additional  approximations  were  needed 
to  successfully  use  the  method.  First,  it  was  found  that  the  computa¬ 
tion  of  the  Bhock  position  and  slope  led  to  instabilities  in  the 
solution,  and  hence  an  approximation  was  introduced  for  the  shock 

slope.  Second,  the  inclusion  of  the  terms  puv  and  Pw  led  to 

s  n 

instabilities  in  the  computacion  near  the  axis  of  symmetry.  These  terms 

entered  into  the  computation  of  the  pressure,  and  the  pressure  influenced 

the  ether  variables  through  the  term  of  the  tangential  momentum 

equation.  To  remove  the  instability,  Davis  and  Chyu  omitted  the 

contributions  of  puv  and  Pw  to  the  tangential  pressure  gradient, 

s  n  n 

dp  dp  °Pm 

That  is,  ^  was  approximated  by  a  function  where  p^  is 

obtained  from  a  simplified  form  of  the  normal  momentum  equation: 


dpT 


epKu 


l  +  e<n 


(*.l) 


The  function  p^  is  identical  to  the  pressure  which  is  computed  in 
the  thin-shock-layer  theory.  Hence,  with  these  two  approximations-- 
on  the  shock  angle  and  on  the  pressure  gradient --the  problem  becomes 
similar  to  the  thin- shock- layer  problem. 

The  specific  procedure  used  by  Davis  and  Chyu  is  outlined  below. 
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1)  The  tangential  momentum  equation  is  programmed  into  the 
implicit  finite-difference  scheme,  and  under  suitable  assumptions  about 
the  shock  slope  and  pressure  gradient,  the  tangential  velocity,  u, 

is  computed  from  the  scheme. 

2)  The  continuity  equation  is  integrated  term  by  term,  and  the 
normal  velocity,  v,  is  then  determined  by  a  simple  numerical  evaluation 
of  an  integral. 

3)  The  pressure,  p,  is  obtained  by  a  similar  integration  of  the 
normal  momentum  equation. 

The  first  tvo  steps  of  this  procedure  are  identical  to  the  procedure 
vhich  has  been  successfully  used  for  solving  the  boundary- layer 
equations  [. l4,  12,  13]  (for  a  variable  density,  the  temperature,  T, 
is  computed  along  vith  u  in  the  first  step).  In  boundary- layer  theory, 
of  course,  the  pressure  is  a  known  function.  Hence  the  addition  of  the 
third  step  is  the  most  obvious  means  of  extending  the  method  to  the 
shock-layer  computations.  The  further  extension  of  this  method  of 
solution  to  flows  of  variable  density  is  considered  in  section  C  and 
is  referred  to  as  method  I. 

It  was  noted  in  Chapter  II  that  Weinbaum  and  Garvine,  in  an 
Investigation  of  the  flow  In  a  laminar  wake,  have  solved  a  set  of 
equations  thatare  identical  to  equations  (2.2)-(2.7)*  The  method  of 
solution  that  was  proposed  In  reference  4o  and  considered  further  in 
reference  39  is  discussed  with  regard  to  its  application  to  blunt-body 
flows  In  section  D  (method  II).  This  method  makes  use  of  the  charac¬ 
teristic  surfaces  of  the  equations  to  evaluate  the  pressure  and  normal 


veloc Ity. 


One  additional  variation  of  the  basic  method  of  implicit  finite- 
differences  is  considered  in  section  E.  In  this  method,  all  four 
variables,  u,  v,  T,  and  p  (or  p),  are  programmed  into  the  implicit 
finite -difference  scheme  (method  III). 

In  the  forms  that  are  considered  in  the  following  sections,  none 
of  these  methods  has  been  found  to  be  entirely  satisfactory.  The  use 
of  the  first  and  third  methods  led  to  instabilities  in  the  computations, 
while  the  second  can  be  seen  to  be  somewhat  inconvenient  for  the  blunt- 
body  problem.  It  should  be  emphasized;  however,  that  only  the  first  of 
these  methods  has  beer  explored  in  detail,  and  the  lack  of  success  in 
applying  these  methods  cannot  be  interpreted  as  conclusive  evidence 
that  such  methods  are  not  feasible.  Hence  the  methods  are  described 
below  in  order  to  indicate  the  type  of  problems  that  are  encountered 
and,  wherever  possible,  to  indicate  the  causes  and  possible  corrections 
for  instabilities  and  errors.  Several  suggestions  have  been  made  for 
the  incorporation  of  certain  features  of  the  first  two  methods  into  the 
method  which  is  described  in  section  E. 

B.  3s sic  Equations  of  the  Implicit  Finite-Difference  Method 

Since  the  basis  of  each  of  the  methods  outlined  above  is  an  implicit 
f inite-diffe *ence  scheme,  we  first  describe  the  basic  principles  of  such 
a  scheme.  All  necessary  equations  have  been  given,  but  seme  details  of 
the  method  have  not  been  included.  For  a  more  detailed  discussion  of 
the  method  and  for  discussions  of  the  merits  of  various  numerical  schemes, 
the  reader  is  referred  to  the  inv?stigations  of  Flugge-Lotz  and 
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Blotther  [l4],  Davis  and  Flugge-Lotz  [ 12],  and  Fannelop  and 
Flugge-Lotz  [13]. 

The  region  "between  the  shock  and  the  "body  is  replaced  with  a  grid, 

which  is  depicted  in  figure  4.1.  The  values  of  the  flow  variables  are 

to  be -evaluated  only  at  the  mesh  points,  which  are  denoted  by  the  sub¬ 
scripts  m  and  N.  If  the  grid  spacing  is  constant  (this  is  not 

essential  but  is  assumed  in  this  chapter),  the  coordinates  are  given 
by  s  =  mAs  and  n  =  NAn.  It  is  further  assumed  that  the  values  of 
the  flow  variables  are  known  at  all  mesh  points  for  s  <  mAs  and  that 
we  wish  to  solve  for  the  values  at  all.  mesh  points  on  the  line 
8  =  (m+l)A  s. 


Figure  4.1.  Finite-difference  grid. 


If  the  partial  derivatives  of  equations  (2 -2 )-(2 -5 )  are  replaced  with 

difference  quotients,  the  nonlinear  partial  differential  equations, 

(2.2)-(2.5)j  for  the  functions  u(s,n),  v(s,n),  etc.  are  replaced  with 

algebraic  difference  equations  for  the  values  of  the  functions  evaluated 

at  the  mesh  points.  The  unknown  quantities  in  this  case  are  denoted  by 

u  v  , ,  etc.  for  N  =  0,1,  ...,N  .  N  denotes  the  grid  point 

immediately  inside  the  shock  and  therefore  is  the  greatest  integer  which 

A 

is  less  than  or  equal  to  ^  . 

The  difference  quotients  that  are  used  for  this  reduction  are  given 


dF  5Fm+l,N  ■  4Fm,N  +  Fm-1,N  .  „/A  2V 
&  =  -  "2 As"  - +  °(AS  > 


(4.2a) 


dF  Fm+1,N+1  "  Fm+1,N-1  .  2, 

K  - - * —  +  > 


(4.2b) 


tli  =  JSt hm.  "  2Fm+^»  *  Fm+1, N-l  +  0(An2} 

dr2  An2 


(4 .?c) 


dF|  dG  =  JL_ 
5n  5n  i.. 


2  [|Gm,N+l  '  Gm,N-l)(Fm+l, 


N+l  “  Fm+1 


,N-1 


+  If  -  f 

j  i,H+l  m,N 


-l)lGm+l, 


N+l  m+1,  N-l 


m,  N+l 


“  Fm,N-l)( 


Gm, N+l  "  Cm,N- 


+  0(As2)  +  0(An2) 


(4. 2d) 
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where  each  of  the  derivatives  has  been  evaluated  at  the  point  (m+l,N). 
Thus,  the  s -derivative  is  replaced  with  a  "backward"  difference  quotient, 
and  the  n-derivatives  are  replaced  with  "central"  difference  quotients. 
The  order  of  magnitude  of  the  errors  associated  with  the  use  of  such 
difference  quotients  has  been  shown  in  the  equations.  The  form  of 
equation  (4. 2d)  was  chosen  so  that  the  unknown  variables  at  m+1  appear 
in  a  linear  manner. 

The  partial  differential  equations  are  evaluated  at  the  grid  points 

(m+l,N),  N  -  1,2, ...,N_,  by  use  of  the  difference  quotients  (4.2). 

s 

This  leads  to  a  system  of  simultaneous  algebraic  equations  for  the 

unknown  variables.  This  coupling  of  the  equations  is  a  result  of 

having  used  a  "backward"  difference  formula  in  equation  (4.2a)  and  is 

the  feature  which  gives  the  method  its  "implicit"  character.  Due  to 

du 

the  presence  of  such  terms  as  u  ^7,  the  algebraic  equations  will  be 
nonlinear.  To  facilitate  the  solution  of  the  equations,  they  are 
generally  linearized  by  the  use  of  extrapolation  formulas  such  as 


■e>  —  2V  -  F 

*m+l,N  m,N  m-l,N 


+  0(As2) 


(4.2e) 


The  solution  that  is  obtained  from  the  linear  difference  equations 
may  itself  be  used  to  linearize  the  equations  for  a  computation  of 
an  improved  solution  if  the  errors  from  the  original  linearization 
prove  to  be  excessive.  Such  an  iteration  procedure  can  be  continued 
until  the  difference  between  successive  solutions  is  as  small  as 
desired.  Hence,  the  solution  may  be  considered  to  be  an  iterative 
solution  of  nonlinear  algebraic  equations  even  though  the  emphasis  in 
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this  chapter  is  on  linear  equations .  It  is  generally  desirable  to  avoid 
such  iterations  since  they  can  lead  to  excessive  computation  times.  In 
previous  applications  of  the  method  [l4,  12,  13],  it  has  been  found 
that  such  iterations  are  not  necessary.  However,  E.  Krause  [22]  has 
considered  several  examples  where  iteration  was  either  necessary  or 
was  more  convenient  than  alternatives. 

The  nonlinear  partial  differential  equations  may  now  be  considered 
to  have  been  reduced  to  a  system  of  coupled,  linear  algebraic  equations. 
For  the  method  considered  in  section  C,  the  tangential  momentum  equation 
and  tne  energy  equation  are  reduced  in  this  manner  and  lead  to  difference 
equations  of  the  form 

A1K  Um+1,N+1  +  Um+1,N  +  C1N  Um+1,N-1  *  D1N  Tm+l,».l 

•".'..l/'Wl,.-!"1.  <4-5a> 

and 


A2  u  , .  _TJ,  +  B2  u  .  +  C2  u  , ,  „  ,  +  D2„  T  .TJ, 

N  m+l,N+l  N  m+1,  N  N  m+1,  N-l  N  m+1,  N+l 


*  Tm+1,N  +  Tm+1,N-1  G2N 


(4ob) 


where  N  =  1,2, ...,N  •  The  coefficients  Al,  Bl, ...,G2  are  functions 

£ 

only  of  the  known  values  of  the  flow  variables  at  s  =  mA s  and  (m-l)As 

and  of  the  extrapolated  values  (or  values  obtained  from  the  preceding 

iteration)  at  s  =  (m+l)As.  Hence,  they  can  be  assigned  numerical 

values.  This  system  of  2N  simultaneous  equations  has  a  special 

s 

"tridiagonal"  form  [l4j  and  consequently  can  be  solved  without  recourse 
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to  the  usual  matrix -inversion  procedures  (which  would  be  inefficient  for 
the  fairly  large  values  of  N  that  are  used) . 

D 

For  this  case  of  two  variables,  u  and  T,  the  solution  can  be 
written  out  in  terms  of  the  coefficients  Al,  Bl, . ..,G2  fairly 
easily.  However,  for  the  case  of  four  variables,  which  is  considered 
in  section  E,  the  equations  become  unwieldy  unless  written  in  matrix 
notation.  Hence,  a  matrix  notation  is  adopted  at  this  point  in  order 
to  handle  both  cases  at  once.  The  function  w(n)  is  defined  as  a 
column  vector  of  the  dependent  variables.  In  sections  C  and  D,  it  has 
two  components, 

fUm+l,Nl 


v(N) 


m+1,  N 


and  in  section  E,  it  has  four  components, 


w(N) 


u 


m+l,N 

rm+l,N 

%i+l,N 


m+l,Nj 

Let  the  number  of  components  be  denoted  by  i.  A  complete  solution  of 
the  difference  equations  has  been  obtained  when  w(n)  has  been 
evaluated  at  N  =  0, 1, . . . ,N +1. 

D 

The  linear  algebraic  equations  that  replace  the  partial  differential 
equations  now  have  the  form 


A(N)  w(N-l)  +  B(.  '  w(N)  +  C(N)  v(N+1)  =  D(N)  (4.4) 
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for  N  =  1,2,... ,N  .  For  each  value  of  N,  the  coefficients  A(N), 

8 

B(N),  and  C(n)  are  i  x  i  matrices,  and  the  coefficient  D(N)  is 

a  column  vector  with  i  elements.  It  is  emphasized  again  that  they 

are  functions  that  can  be  evaluated  from  known  data,  from  extrapolated 

data,  or  from  data  obtained  in  the  preceding  iteration.  For  future 

reference,  it  is  noted  that,  in  the  case  of  i  =  2,  one  component  of 

D(N)  is  a  function  of  the  tangential  pressure  gradient,  ^  ,  evaluated 

at  (s,n)  =  ((m+l)A  s,  NAa).  Equation  (4.4)  is  a  linear,  second-order 

difference  equation  and  represents  N  equations  for  the  N  +2  unknowns, 

v(0), ...,w(N  +1).  Two  additional  equations  are  obtained  from  boundary 
8 

conditions.  These  boundary  conditions  may  take  the  form  of  first-order 
difference  equations  without  altering  the  tridiagonal  form  of  the 
equations.  Hence, 


B(C)  w(0)  +  C(0)  w(l)  =  D(0) 


(4.5a) 


and 

A(Nb+1)  w(Ng)  +  B(Nb+1)  w(Ns+1)  =  D(NS+1)  (4.5b) 

are  the  acceptable  forms  of  the  boundary  conditions.  In  practice,  the 

requirement  that  the  boundary  conditions  take  this  form  does  not  pose 

a  restriction.  As  an  example,  consider  a  boundary  condition  of  the 

dTl 

form  =  0.  The  derivative  must  be  replaced  with  a  difference 

quotient.  If  a  two-point  form  is  used,1'  T(l)-T(0)  =  0,  the  boundary 


The  notation  that  was  adopted  to  describe  the  vector  w  is  also 
frequently  used  for  the  other  variables.  In  this  notation,  the  sub¬ 
script  m+1  is  dropped,  and  it  is  understood  that  the  function  is 
evaluated  at  s  =  (m+l)As.  In  addition,  the  variable  is  considered  to 
be  a  function  of  the  index,  N,  and  this  functional  dependence  is  written 
exp lie itly- 
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condition  has  the  form  of  equation  (4.5a).  Unfortunately,  the  two-point 

form  generally  does  not  result  in  a  sufficiently  accurate  description 

of  this  boundary  condition.  Use  of  a  three-point  difference  quotient, 

*>  -  3T(°)  ^T(1)  'l'T(2)  _  o,  is  more  accurate  but  results  in  a 

dn|n=0  2An 

boundary  condition  of  the  form 

B(0)  w(0)  +  C(0)  w( 1 )  +  w(2)  =  D(0)  , 

which  does  not  conform  to  the  desired  tridlagonal  form.  However,  we 
know  from  the  theory  of  linear  algebra  that  we  can  add  to  any  equation 
a  linear  combination  of  the  other  equations  without  altering  the  solu¬ 
tion.  Hence,  we  combine  the  "non -conforming"  boundary  condition  with 
equation  (4.4)  evaluated  at  N  =  1  in  such  a  way  that  w(2)  is 
eliminated.  Thus 

[A( l)  +  C ( l)  B(0)]  w(0)  +  [B ( 1 )  +  C(l)  C(0)]  w(l) 

=  D(l)  +  C (l)  D(0) 

would  contain  the  desired  information  that  5T(0)  -  4t(1) +T(2)  =  0 
and  would  still  conform  to  the  necessary  tridiagonal  form. 

The  solution  to  the  second-order  difference  equations,  (4.4)  and 
(4.5),  is  given  by  a  first-order  difference  equation 

w(N)  =  H(N)  +  K(N)  w(N+1)  ,  N  =  0,1,. ..,N  (4.6) 

(for  proof  of  this  statement,  see  reference  14).  The  coefficients,  H 
(a  column  vector)  and  K  (an  i  X  i  matrix),  are  evaluated  from 
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recursion  formulas,  which  are  obtai:  d  as  follows.  Replace  N  with 
N-l  in  equation  (4.6)  to  obtain 

w(N-l)  =  H(N-l)  +  K(N-l)  w(N)  ,  N  =  1,2, ...,N  +1 

B 

This  equation  can  be  us?d  to  eliminate  w(N-l)  from  (4.4): 

[B(N)  +  A(N)  K(N-l)3  w(N )  +  C(N)  w(M+l)  =  D(N)  -  A(N)  HU-!) 

Solving  this  equation  for  w(n)  and  comparing  the  result  to  equation 
(4.6),  we  obtain 

H(N)  =  [B(N)  +  A(N)  K(N-l)]’1  (D(N)  -  A(N)  H(K-l)3  (4. 7a) 

and 

K(N)  =  -lB(H)  +  A(N)  K(N-l)]"1  C(N)  (4.7b) 

The  values  of  H  and  K  at.  the  surface  of  the  body  can  be  obtained  by 
inspection  when  the  equation  w(0)  =  H(0)  +  K(0)  w(l)  is  compared  to 
the  boundary  condition  (4.5a),  i.e., 

H(0)  =  B(O)"1  D(0) 

and  (4.8) 

K(0)  =  -B(O)"1  C(0) 

Recursive  application  of  equations  (4.7)  for  N  =  1,2, ...,N  yields 

s 

the  values  of  H  and  K  at  all  grid  points  across  the  shock  layer. 

Algebraic  combination  of  the  equation  w(N  )  =  H(N  )  +  K(N  )  w(N  +l) 

with  the  boundary  condition  (4.5b)  yields  the  value  of  w(N  +l)-  Then 

& 

a  recursive  application  of  equation  (4.6)  for  N  =  N  ,  N  -1, — ,1,0 

&  & 

yields  the  desired  solution  to  the  system  of  equations.  It  should  be 
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noted  that  a  straightforward  solution  to  equations  (4.4)  and  (4.5)  by 
standard  matrix  inversion  techniques  would  require,  at  each  s-step, 
the  inversion  of  a  square  matrix  of  dimensions  (n  +l)i  X  (N  +l)i. 

8  S 

This  would  require  prohibitive  computation  times  if  N  is  large. 

5 

The  solution  given  by  equations  (4.6)  and  (4.7),  on  the  other  hand, 

requires,  at  each  s-step,  N  inversions  of  an  i  x  i  matrix.  Since 

s 

i  is  never  greater  than  four  for  the  shock  layer  equations,  these 

operations  can  be  efficiently  carried  out  for  large  values  of  N  . 

The  use  of  matrix  notation  has  resulted  in  a  rather  concise 

description  of  the  solution  in  equations  (4.6)  and  (4-7).  The  matrix 

fomulas  can  be  used  intact  in  a  computer  program  if  one  makes  use  of 

arrays  (or  subscripted  variables)  in  the  machine.  However,  this 

generally  results  in  inefficient  operation  of  the  program  due  to  an 

increase  in  access  time  to  the  stored  data.  It  should  be  noted  that 

it  is  necessary  to  store  the  values  of  w(n),  H(n),  and  K(N)  at  each 

grid  point  across  the  shock  layer  (A(n),  B(n),  C(n),  and  D(N)  are 

needed  only  once  and  hence  need  not  be  stored).  This  would  result  in 

the  function  K(N)  being  a  three-dimensional  array  in  the  machine-- 

with  N  +1  elements  in  one  dimension  and  i  elements  in  the  other  two. 
s 

The  access  time  to  data  that  is  stored  in  such  higher-dimensional  arrays 
is  usually  significantly  larger  than  that  for  one-dimensional  arrays. 

For  the  case  when  only  u  and  T  are  programmed  into  the  implicit 
finite-difference  scheme,  i  =  2,  and  the  solutions  given  by  equations 
(4.6)  and  (4.7)  can  easily  be  written  out  explicitly  in  terms  of  the 
coefficients  of  equation  (4.2)  and  the  individual  elements  of  H  and  K. 
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These  expressions,  though  less  concise  than  those  given  above,  are  more 
efficient  for  the  actual  computation  of  the  solution.  The  reader  is 
referred  to  reference  14,  12,  or  for  these  expressions.  To  a  certain 
extent,  this  can  be  done  for  the  case  when  i  =  4,  but  the  expressions 


l 


are  lengthy,  and  unless  the  computation  time  is  a  critical  factor,  the 
extra  effort  needed  to  obtain  the  expressions  may  not  be  considered 
worthwhile. 

Using  these  basic  expressions  for  the  solution  of  the  implicit 
finite-difference  equations,  we  can  now  consider  specific  methods. 


C .  Method  I 

The  computation  procedure  that  has  been  most  thoroughly  examined 
is  one  that  is  similar  to  the  method  used  by  Davis  and  Chyu  in  their 
study  of  the  constant-density  shock  layer.  In  this  method,  only  the 
tangential  momentum  equation  and  the  energy  equation  have  been  reduced 
to  linear  difference  equations  of  the  form  (4.5).  The  specific  expres¬ 
sions  for  Al, Bl, ...,G2  have  been  listed  in  Appendix  A. 

The  boundary  conditions  which  complete  the  set  of  equations  are 
obtained  from  the  slip  and  temperature -jump  conditions  at  the  wall, 
equations  (2.10)  and  (2.11),  and  from  the  Rankine-Hugoniot  conditions 

for  u  and  T,  equations  (2.15)  and  (2.17).  In  both  (2.10)  and  (2.11), 

dui  5u(0)-4u(l)  +  u(2)  _  ,  ,  ,  , 

we  set  ‘ — 2An  - nas  alrea<iy  13  -en 

explained  how  the  resulting  three-point  boundary  condition  can  be 

reduced  to  the  necessary  two-point  form.  The  function 

T(0)  ,  which  appears  in  both  (2-10)  and  (2.1l),  must 


p(o) 


Jf 
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be  evaluated  by  extrapolation  (equation  (4.2e))or  from  a  previous 
iteration  in  order  to  have  a  linear  boundary  condition.  The  shock 
conditions  are  somewhat  more  difficult  to  impose  since  the  shock  loca¬ 
tion  will  not,  in  general,  coincide  with  a  grid  point.  Hence,  it  is 
necessary  to  use  interpolation  formulas  in  order  to  apply  the  Rankine- 

Hugoniot  conditions.  Remembering  that  N  is  the  largest  integer 

_  s 

which  is  less  than  or  equal  to  ^  ,  we  define  the  shock  location  as 

A=(Ns+Qs)An  (4.9) 

Hence,  0  <  a  <1.  The  values  of  the  dependent  variables  behind  the 

“  5 

shock,  which  we  denote  by  wA(0),  are  given  as  functions  of  0  by 
equations  (2.13)  and  (2.17).  Suitable  interpolation  formulas  that  can 
be  used  as  boundary  conditions  in  the  difference  scheme  are 


Two -point  or  linear  interpolation 
w(Ns)(l-OJ  +  v(Ns+l)ag  =  wA(c)  (4.10a) 

and 

Three-point,  or  quadratic  interpolation 
a 

w(Ns+l)  —■  (HO.)  +  w(Ng)(l-CT) 
a 

-  w(Ns-l)  -f  (1-0.)  =  wA(0)  (4.10b) 


The  latter  formula  requires  a  reduction  to  the  proper  two-point  form. 
This  can  lead  to  some  rather  unpleasant  algebra,  and  therefore,  in  one 
case  to  be  considered  later,  the  simpler  equation  (4.10a)  was  used. 
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For  the  application  of  these  boundary  conditions  and  thus  for  the 
computation  of  u  and  T  from  the  difference  equations,  it  is 
necessary  to  know  the  values  of  A  and  0  .  This  can  be  handled  in 
several  ways.  For  the  moment,  we  treat  A  and  0  in  the  same  manner 
that  the  flow  variables  are  handled.  That  is,  we  extrapolate  their 
values  from  previous  data,  and  using  this  extrapolated  data,  we  can 
compute  new  values  of  A  and  0  at  a  later  point  in  the  procedure. 
Actually,  only  one  of  the  variables,  A  or  0  ,  needs  to  be  extrapolated* 

the  remaining  quantity  can  be  obtained  from  a  difference  representation 
of  equation  (2.19).  It  will  later  be  shown  that  it  is  possible  to 
compute  A  simultaneously  with  the  computation  of  u  and  T  without 
using  this  extrapolation.  However,  for  the  moment,  we  consider  u  and 
T  to  be  computed  from  equations  that  make  use  of  approximate  values 
of  A  and  0 

The  normal  velocity,  v,  is  obtained  from  the  continuity  equation 
after  the  computation  of  u  and  T  has  been  accomplished.  An  Integra¬ 
tion  of  the  terms  of  equation  (2.2)  yields 

1  C  ^  1 

p(s,n)v(s,n)  =  - i -  I  ((r+ensin  0)JPu]  dn 

(r+ensin  0)*^(l  +  K€n)  Jo  8 

(4.11) 

In  order  to  evaluate  v  from  this  equation,  it  is  necessary  to  use  an 
extrapolated  value  for  p.  This  is,  however,  consistent  with  the 
linearization  of  the  difference  equations  for  u  and  T.  The  boundary 
condition  (2.12)  has  already  been  incorporated  into  (4.11).  The 


115 


numerical  evaluation  of  the  integral  in  equation  (4.1l)  is  carried 
out  by  the  use  of  a  formula  similar  to  Simpson's  rule1: 


J(H-l)An 


F  dn  =  —  [5F(N-1)  +  8F(N)  -F{N+l)] 


(4.12) 


This  formula  is  derived  by  replacing  F(n)  with  a  quadratic  curve 
which  passes  through  the  three  points,  F(N-l),  F(n),  and  F(N+l). 

At  this  point  in  the  computation,  an  improved  value  for  the  shock 
position.  A,  may  be  obtained  simply  by  observing  where  the  integrated 
value  of  v  reaches  the  Rankine-Hugoniot  value  given  by  equation  (2.l4). 
One  may,  of  course,  now  repeat  the  computation  of  u  and  T  using 
the  new  value  of  A  (and  hence  0  )  until  the  value  of  A  is  as 
accurate  as  desired.  Alternatively,  one  can  proceed  to  the  computation 
of  the  pressure  and  then  iterate  on  all  of  the  approximated  variables 
simultaneously. 

Once  the  values  of  u,  T,  and  v  have  been  obtained  at  each  grid 
point,  the  pressure  can  be  evaluated  from  an  integration  of  the  normal 
momentum  equation.  If  the  extrapolated  value  of  p  is  used,  the 
pressure  can  be  evaluated  from 

r n  Teuv  2 

p(s,n)  -  p(0)  -  a  ^  +  €Wn  -  (4.13a) 


Although  this  is  consistent  with  the  linearization  of  the  finite-difference 


Simpson's  rule  is  not  used  because  it  computes  the  integral  from 
(N-l)An  to  (N+l)An  and  therefore  would  provide  the  value  of  the 
integral  only  at  alternate  grid  points  across  the  shock  layer. 
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scheme,  the  approximate  value  of  p  need  not  be  used  since  the  state 
equation  can  be  used  to  eliminate  P  from  equation  (2.-.^.  This  leads 
to 


p(s,n)  =  p(0)  exp 


['f p(n)  H 


(4.13b) 


where 


P(n) 


_X_  £ 

Y-l  T 


r  euv. 


,  ,  J  +  €W 

|  1+K'£n  n 


k£  1 
l-Hcen  J 


Since  all  the  variables  that  appear  in  the  expression  for  P  have 
been  computed,  equation  (4.13b)  can  be  evaluated  without  use  of  the 
extrapolated  data.  In  both  (4.13a)  and  (4.13b),  the  value  of  p(0) 
is  chosen  so  that  the  pressure  satisfies  the  Rankine-Hugoniot  condition, 
equation  (2.15). 

Several  possible  variations  in  the  details  of  the  computation 
procedure  have  already  been  noted.  There  are  other  variations  that 
have  been  investigated  but  are  not  described  here.  In  general,  when 
both  the  shock  shape  and  the  tangential  pressure  gradient,  ^  ,  were 
computed  at  each  s-step,  the  computed  values  become  meaningless  within 
a  few  steps.  In  other  cases,  the  procedure  did  not  exhibit  such  strong 
instabilities  but  yielded  incorrect  results  nevertheless.  As  noted 
above,  such  difficulties  were  not  entirely  unexpected  and  can  be 
attributed  to  the  determination  of  the  shock  shape  and  to  the  occurrence 
of  uvfi  and  wn  in  the  normal  momentum  equation.  Qualitatively,  the 
difficulty  with  the  normal  momentum  equation  is  as  follows.  Equations 
such  as  (4.2a)  or  (4.2b)  are  used  to  numerically  evaluate  certain 


117 
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derivatives.  When  an  erro^  exists  in  thu  value  of  the  function  T.  the 

use  of  the  difference  quotients  results  in  a  magnification  of  the  error 

in  the  derivative  of  P  (due  to  the  small  value  of  As  or  An  that 

is  used  in  the  denominator).  In  the  basic  implicit  finite -difference 

scheme,  this  does  not  occur  since  equations  (4.2)  are  only  used  formally 

but  are  not  numerically  evaluated.  However,  in  the  procedure  described 

above,  ^  is  evaluated  from  (4.2a)  and  is  used  in  the  determination 

of  u  and  T.  The  resulting  value  of  u  is  used  in  (4.2a)  to 

du 

evaluate  ^7  ,  which  is  in  turn  used  in  equation  (4.1l)  to  evaluate 
v.  From  v,  we  determine  both  ^  and  ^  ,  and  these  values  are 
used  in  equation  (4.13)  to  determine  p.  This  value  of  p  then  goes 
back  into  the  finite-difference  scheme  as  At  each  step  the  evalua¬ 

tion  of  the  derivatives  results  in  a  magnification  cf  the  errors,  and 
these  errors  are  eventually  fed  back  into  the  scheme.  In  the  applica¬ 
tion  to  boundary-layer  flows,  such  a  "feedback"  does  not  occur  since 
the  pressure  is  a  known,  specified  function.  Even  in  the  thin-shock- 
layer  flow, where  the  pressure  is  computed,  the  feedback  does  not  occur 
since  uvg  and  wn  do  not  appear  in  the  normal  momentum  equation. 

We  thus  consider  modifying  the  equations  in  the  manner  of  Davis 
and  Chyu  [J.l] .  This  modification  can  be  considered  to  be  a  first  step 
to  the  more  general  analysis;  the  numerical  scheme  must  be  able  to 
handle  the  simpler  equations  before  they  can  be  extended  to  the  general 
case.  Thus  the  shock  angle,  0  ,  is  specified  as  a  known  function  for 
the  evaluation  of  the  Rankine-Hugoniot  conditions.  Haring  the  course 
of  the  computation,  the  shock  position,  A,  is  calculated.  The 
geometric  relation  between  A  and  0  ,  equation  (2.19),  cannot  now 


118 


be  applied  but  can  only  be  used  a  posteriori  to  determine  whether  the 
results  are  consistent  with  the  assumptions.  In  reference  11,  for  the 
constant -density  flow,  0  was  specified  as  being  equal  to  9,  and 
the  resulting  value  of  A  was  found  to  be  consistent  with  this  assump¬ 
tion  (i.e.  A  ~  const).  For  a  variable  density,  the  modification  of 
the  normal  momentum  equation  can  be  done  in  two  ways.  First,  the 

terms  proportional  to  uv  and  w  may  be  eliminated  entirely  from 

8  n 

the  problem,  and  the  pressure  p  is  approximated  by  the  thin-shock- 

layer  pressure  p^.  The  problem  Is  then  reduced  to  one  that  represents 

a  slight  extension  of  the  investigation  of  H.K.  Cheng.  The  additional 

terms  of  0(c)  that  are  included  in  equations  (2.3)-(2.5)  should  not 

cause  any  major  difficulty  in  the  computation.  Despite  the  fact  that 

the  Bhock  angle  is  now  assumed  to  be  known  (and  therefore  the  boundary 

values  at  the  shock  are  known),  the  treatment  of  the  shock  position 

is  still  somewhat  more  general  than  that  used  by  Cheng  since  the  mass 

flow  in  the  shock  layer  is  still  not  s  specified  quantity.  The  omission 

of  the  terms  puv  and  pw  results  in  a  substantial  reduction  of 
s  n 

the  pressure,  and  therefore  of  the  density,  near  the  body  surface. 

This  reduction  in  the  density  alters  the  flow  field  considerably  and, 
in  an  extreme  case,  can  lead  to  the  occurrence  of  the  zero-pressure 
point  which  is  common  in  thin-shock-layer  theory.  Hence,  a  second 
possibility  that  Is  analogous  to  the  analysis  made  by  Davis  and  Chyu 
is  also  considered.  The  complete  pressure  is  used  in  the  computation 
of  the  density,  and  the  pressure  gradient  ^  is  approximated  by 
dpT 

ST  >  where  PT  iS  the  thin -shock -layer  pressure  of  equation  (4.1). 
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This  should,  result  in  a  more  realistic  description  of  the  density  field, 
but  the  influence  of  puvs  and  pwn,  which  is  fed  back  into  the 
finite -difference  scheme  through  the  density,  may  be  large  enough  to 
cause  instabilities. 

The  first  of  these  possibilities  is  now  considered.  The  computa¬ 
tions  proceed  as  follows.  The  solution  of  u,  T,  and  v  is  based 
on  extrapolated  data  as  described  earlier.  The  shock  location  is 
determined  by  a  requirement  that  the  integrated  value  of  v  be  eaual 
to  the  Rankine-Hugoniot  value.  Since  0  has  been  taken  to  be  a  known 
value,  this  step  is  now  particularly  simple;  and  since  equation  (2.19) 
no  longer  acts  as  a  constraint,  the  computed  value  of  A  may  undergo 
rather  rapid  changes  in  value.  Using  the  newly  computed  values  of  u, 
v,  T,  and  A,  the  pressure  is  obtained  from  equation  (4.13b)  (with 
the  terms  proportional  to  uv  and  w  omitted,  of  course).  As 
many  iterations  as  desired  may  be  obtained  at  each  step. 

This  procedure  has  been  applied  to  the  flow  problem  that  was 
considered  in  Chapter  III,  section  D.l,  i.e.,  y  =  1.4,  M  =  10. 

Res  ~  100'’  a  =  b  and  H  =  T1^2.  These  computations  were 

started  from  initial  data  (at  s  =  0.0  and  s  =  0.02)  that  were  obtained 
from  the  first  truncation  problem  of  Chapter  III  (the  second  truncation 
problem  had  not  been  solved  when  these  computations  were  made).  Thus 
the  shock  angle  0  was  set  equal  to  the  body  angle  9.  It  has 
been  found  that  the  computed  shock  description  is  a  sensitive  indicator 
of  the  stability  of  the  computations.  Thus,  the  computed  value  of  A 
is  shown  in  figure  4.2.  The  discontinuous  jump  in  value  that  occurs 
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NOIlVOCn  XOOHS  ‘2 


Shock  standoff  distance  computed  from  method  I  using  the  thin-shock-layer  approxi 
raation  for  the  tangential  pressure  gradient  and  the  density:  y  ■  1.4,  M  *  10, 

Re  =  100,  b  =  0.6,  o  =  0.7,  and  =  1/2. 


at  the  first  computation  step  is  of  particular  interest.  It  was  found 
after  these  computations  were  made  that  the  initial  data  were  incorrect. 

The  correct  first-truncation  value  of  A  has  also  been  shown  in  figure 
4.?.  It  can  be  seen  that  the  finite -difference  scheme  yields  the  correct 
value  of  A  even  when  started  from  an  incorrect  value.  As  s  increases, 
the  number  of  iterations  that  are  required  to  obtain  "convergence"  of 
the  pressure  increases.  In  addition,  the  computed  values  of  A  become 
less  consistent  with  the  assumption  that  0  =  0>  at  s  =  0.34,  the 
application  of  equation  (2.19)  yields  0-0  =  11°.  Thus  there  is  little 
point  in  continuing  the  calculations  to  larger  values  of  6.  This 
inconsistency  in  the  computation  would  no  doubt  be  removed  by  the 
specification  of  the  mass  flow  in  the  shock  layer  (as  is  done  in  the 
thin-shock-layer  computation).  If  this  were  done,  the  resulting 
problem  would  represent  an  extension  of  the  results  of  H*K.  Cheng 
to  a  more  uniform  treatment  of  the  second-order  boundary-layer  effects. 

It  should  be  noted  that  if  the  shock  angle  is  allowed  to  vary  in 
the  shock  conditions  and  if  equation  (2.19)  is  used  to  relate  A  and  0 
in  the  above  calculations,  the  computed  value  of  A  very  slowly  decreases. 
A  similar  result  had  been  obtained  earlier  when  considering  the  complete 
normal  momentum  equation.  This  is  typical  of  the  difficulties  associated 
with  attempting  to  compute  the  shock  shape. 

We  now  consider  the  second  treatment  of  the  normal  momentum 

equation.  The  calculations  include  the  effect  of  the  terms  uv  and 

s 

vvr  in  the  computation  of  the  density  in  order  to  improve  the  descrip¬ 
tion  of  the  density  field.  In  addition,  these  computations  use  an 
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improved  means  of  calculating  the  shock  position.  Instead  of  using 
the  iterative  method  described  earlier,  the  shock  position  is  obtained 
simultaneously  with  the  solution  of  the  finite-difference  equations  for 
u  and  T  as  described  in  the  follovinp  paragraph- 

Earlier  in  this  chapter,  it  was  stated  that  the  calculation  of 
v(N  +1)  was  achieved  by  combining  the  equation  w(N  )  =  fl(N  )+K(N  )w(N  +l) 
with  the  boundary  condition  (4.5b).  Later,  it  was  observed  that  inter¬ 
polation  formulas  (equations  (4.10))  were  used  to  obtain  the  boundary 
condition  (4.5b).  This  led  to  the  appearance  of  A  and  0  in  the 
coefficients  of  equation  (4.5b).  This  was  previously  handled  by  using 
extrapolated  data  for  the  shock  parameters.  However,  if  one  additional 
equation  could  be  obtained  that  related  w(N  ),  w(N  *-l),  and  A,  then 

S  6 

A  could  be  computed  simultaneously  with  w(N  +1),  and  the  necessity 

s 

of  iterating  on  the  shock  position  would  be  eliminated-  Such  an 
equation  can  be  obtained  from  the  integral  form  of  mass  conservation, 
equation  (2.l3).  Let 

<  r  NAn 

l(N)  =  2°e  I  Pu(r+ensin  0)J  dn  (4.l4a) 

Jo 

It  is  shown  in  Appendix  B, for  the  case  of  1  =  2,  that 

l(N)  =  Ar(N)  +  Bj.(N)u(n)  +  CI(N)T(N)  ( 4 .  l4b ) 

where  A^.(n),  Bj(n),  and  Cj.(n)  are  scalar  quantities  that  can  be 
numerically  evaluated  along  with  H(N)  and  K(n);  i.e.,  A^N),  B^N), 

and  Cj.(n)  can  be  evaluated  from  recursive  relations  starting  at 
N  =  0.  Thus  equation  (2.18)  can  be  written  as 
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(r  +  eAn(N6+a>in  0)J+1  =  (l-a^A^ )  + 

BJ(Ns)  w(Ns))  +  as(AI(Ns+l)  +  BI(Ns+l)  w(Ns+l))  (4.15) 

where  the  linear  interpolation  formula  (equation  (4.10a))  has  been  used 
to  evaluate  the  integral  at  the  shock.  Equation  (4.15)  is  the  additional 
equation  that  permits  the  solution  of  A.  The  three  equations  are  non¬ 
linear  and  thus  must  be  solved  by  some  appropriate  numerical  technique. 
Once  the  values  of  A,  w(N  ),  and  w(N  +l)  have  been  obtained,  the 
solution  proceeds  as  described  before. 

With  the  computation  of  A  handled  in  this  manner,  the  Rankine- 
Hugoniot  condition  on  the  normal  velocity,  v,  should  become  superfluous. 
As  noted  in  the  preceding  chapter,  equation  (2.18)  is  not  an  independent 
condition;  it  can  be  used  only  to  replace  one  of  the  boundary  conditions 
which  describe  the  mass  flux  across  the  boundary  Hence  the  integrated 
value  of  v  (equation  (4.11))  should  now  automatically  satisfy  the 
Rankine-Hugoniot  condition.  In  practice,  of  course,  the  linearization 
of  the  finite -difference  equations  and  the  use  of  an  extrapolated  value 
of  p  in  equations  (4.1l)  and  (4.l4a)  prevent  this  from  being  exact  -- 
there  will  always  be  a  small  difference  between  the  integrated  and 
the  Rankine-Hugoniot  values  of  v. 

After  the  pressure  is  determined  from  equation  (4.13b),  iterations 
on  the  values  of  all  the  variables  can  be  obtained.  It  was  found  that 
for  these  calculations,  the  iterations  did  not  significantly  alter  the 
results,  but  did  tend  to  smooth  the  results  somewhat  near  the  axis. 
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The  computed  value  of  A  is  again  used  as  an  indication  of  the 
stability  of  the  scheme.  However,  it  is  convenient  to  make  use  of 
equation  (2.19)  and  to  present  the  results  in  the  form  of  the  shock 
angle.  It  should  be  noted  that  a  small  oscillation  in  the  value  of 
A  results  in  a  rather  substantial  oscillation  in  the  value  of  0 . 

Figure  4.3  illustrates  the  results  for  a  variety  of  assumptions. 

l)  First,  consider  the  cases  for  which  the  computations  were 
started  from  initial  data  obtained  from  the  first -truncation  results. 

As  before,  It  was  found  that  the  computation  of  the  shock  angle  lead 
to  instabilities.  Hence,  the  results  shown  in  figure  4.3  are  based 
on  an  assumed  value  of  0=0  (the  solution  of  the  equations  provides  a 
computed  value  of  0  which  is  the  value  shown  in  figure  4. 3). 

,  dPT 

a)  When  the  pressure  gradient  is  computed  at  each  step, 

instabilities  appear  Immediately.  This  is  illustrated  by  one 

curve  in  figure  4.3.  Thus  the  procedure  is  less  stable  than  it 

was  when  the  influence  of  the  terms  uv  and  vv  was 

s  n 

eliminated  entirely.  This  must  be  a  consequence  of  the  inclu¬ 
sion  of  the  influence  of  uvg  and  wn  in  the  computation  of 
the  density. 

b)  The  strength  of  this  influence  is  evaluated  by  considering 

the  pressure  gradient  to  be  a  specified  function  (equal  to  the 

series -truncation  results).  Hence  the  only  source  of  instability 

due  to  the  terms  uv„  and  w  is  in  the  computation  of  the 

s  n 

density.  The  calculations  now  proceed  to  about  s  =  0.2  before 
instabilities  cause  the  computation  to  break  down  (see  figure 

4-3). 
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FINITE  DIFFERENCE  RESULTS. 


Figure  4.3 

Shock  angle  computed  from  method  I  under  the  approximations  explained  in  paragraphs  ia 


c)  The  onset  of  these  instabilities  can  be  delayed  by  averaging 
the  computed  and  the  extrapolated  values  of  the  pressure: 

P  -  (l  h)  Pgxtrapol^ecL  +  ^  ^calculated  ^ 

At  s  =  0,  h  is  taken  to  be  unity.  The  computations  proceed 
for  increasing  s  until  they  become  unstable.  The  value  of  h 
is  then  decreased,  and  the  calculations  are  started  at  a  point 
just  before  the  onset  of  the  instability.  By  repeating  this 
process  whenever  necessary,  the  calculations  have  been  carried 
to  s  =  0-50,  and  this  result  is  also  shown  in  figure  4.3.  The 
final  value  of  h  is  0.10.  The  instabilities  which  appear  at 
intervals  along  the  curve  are  not  shown  but  ire  quite  similar 
to  the  one  that  has  been  shown  at  s  =  0.2  in  the  figure. 

It  would  seem  that  at  best  these  numerical  computations  could  only 
reproduce  the  truncation  results  since  a  large  number  of  constraints 
have  been  placed  on  the  numerical  procedure.  However,  the  computed 
values  of  0  closely  match  the  second-truncation  result  for  0  even 
though  the  computed  values  are  based  on  initial  data  and  a  pressure 
gradient  that  are  obtained  from  the  first  truncation  and  the  first- 
truncation  assumption  that  0=9  has  been  used  in  the  evaluation 
of  the  shock  conditions.  Hence,  m  this  one  respect,  the  finite- 
difference  scheme  is  yielding  "new" information. 

2)  Also  shown  in  figure  4.3  are  the  results  of  using  the  second- 
truncation  computations  for  initial  data,  for  the  pressure  gradient, 
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and  for  the  shock  angle.  These  computations  have  been  carried  out  to 
s  =  0.62  by  averaging  the  pressure  as  described  in  paragraph  (ic) 
above.  Again,  the  value  of  h  was  gradually  reduced  to  a  final  value 
of  0.10.  No  attempt  has  been  made  to  determine  how  far  the  calcula¬ 
tions  can  be  carried  in  this  manner.  At  the  largest  value  of  s,  most 
of  the  "inviscid",  outer  flow  has  become  supersonic.  This  is  illustrated 
in  figure  4.4  where  the  shock  and  the  sonic  line  have  been  depicted. 

Shown  for  comparison  is  the  sonic  point  on  the  body  as  computed  from  an 
inviscid  theory  (this  value  was  obtained  from  reference  12  and  is  the 
result  of  an  analysis  by  H.  Lomax  of  NASA).  Due  to  the  presence  of  the 
viscous  layer,  the  sonic  line  will  turn  and  roughly  parallel  the  body  as 
s  increases  further.  The  extent  of  the  viscous  layer  has  been  illus¬ 
trated  by  showing  the  tangential  velocity  profile  at  s  =  0.4. 

Returning  to  a  consideration  of  figure  4.3,  we  can  see  that  the 
computed  value  of  0  deviates  only  slightly  from  the  assumed,  second- 
truncation  value.  It  was  previously  noted  that  the  finite -difference 
result  for  0  was  obviously  more  accurate  than  the  first-truncation 
value  or  which  the  computations  were  based.  Whether  this  is  true  of 
the  present  results  can  be  determined  in  the  following  manner. 

In  Chapter  III,  it  was  seen  that  an  inability  to  adequately  deter¬ 
mine  the  shock  shape  was  the  major  cause  of  inaccuracy  in  the  series- 
truncation  analysis .  We  now  use  the  finite-difference  result  found 
here  to  determine  the  value  of  0^,  which  in  Chapter  III  was  arbitrarily 
set  to  zero.  That  is,  we  determine  0^  such  that  0^sin  3  +  0^  sin^s 
(see  equation  (3 -3b))  provides  a  least-squares  fit  to  the  finite -difference 
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O  FIRST-ORDER  HWISCID 
THEORY 


The  shock  and  the  sonic  line  computed  from  method  I 
Y  =  1.4,  Mw  =  10,  Reg  =  100,  b  =  0.6,  or  =  0-7,  and 


result  for  0-  0  .  We  are  not  free  to  choose  the  value  of  0^  as  it  is 

entirely  dictated  by  the  solution  of  the  second -truncation  problem  as 

described,  in  Chapter  III.  Further,  the  values  of  0^  and  0^  are 

mutually  dependent:  0^  sin  s  increases  above  the  value  shown  in 

figure  4.3  as  the  value  of  0^  is  increased.  Hence  the  leaet-squares 

3 

fit  of  01  sin  s  +  0^  sin  s  is  not  as  accurate  as  might  be  expected: 

0-0  =  0.152  sin  s  +  0.112  sin^ s,  and  the  root-mean-square  error  is 
approximately  0.006  compared  with  0-0  =  0.111  at  s  =  0.6.  However,  the 
shear  stress  result  is  now 

t/e  -•  1.23  sin  s  -  0-74  sin^s  , 

which  is  considerably  better  than  the  previous  second-truncation  result 
given  by  equation  (3-24b)  and  illustrated  in  figure  3*1»  It  is  certainly 
not  recommended  that  the  finite-difference  solutions  be  used  to  improve 
the  series  truncation  solutions  in  this  manner;  but  the  above  result 
does  show  that  the  finite-difference  solution  for  0  is  more  accurate 
than  the  second-truncation  result  on  which  the  finite -difference  computa¬ 
tions  were  based. 

Despite  this  ->ne  positive  result,  it  is  quite  clear  that  the  com¬ 
pressible-flow  computations  are  far  more  unstable  than  were  the 
constant -density  computations  of  reference  11.  In  fact,  it  appears  that 
this  method  of  computation  is  not  suitable  for  the  computation  of  the 
flow  in  the  shock  layer;  a  better  means  of  handling  the  computation  of 
the  normal  velocity  and  the  pressure  is  needed. 
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T>.  Method  II 


4 


One  such  possible  method  for  an  improved  calculation  of  the  pressure 
and  normal  velocity  makes  use  of  the  method  of  characteristics.  This 
method  has  several  advantages  over  the  one  described  previcusly.  It 
was  developed  by  Weinbaum  and  Garvine  [4o]  for  use  in  the  laminar  wake 
flow.  Its  application  to  the  blunt-body  problem  is  now  considered. 

In  Chapter  II  it  was  noted  that  there  are  two  "non -parabolic" 
characteristic  surfaces  of  the  partial  differential  equations  (2.2)-(2.5). 
These  surfaces  are  described  by  equation  (2.9): 


-i—  X  +  /5  x 

l+€<n  ds  u  ~  eu  1,2 


(4.17) 


These  two  characteristics  can  be  traced  to  the  continuity  equation  and 
the  normal  momentum  equation,  and  they  are  the  characteristics  which  are 
associated  with  the  determination  of  the  derivatives  of  pressure  and 
normal  velocity.  Hence  we  now  consider  equations  (2.3)  and  (2.5). 

First,  introduce  two  new  coordinates  £  -  |{s,n)  and  T|  =  T|(s,n). 
In  general, 

3  _  3  .  Stj  3 

Ss  5s  5F  Ss 

and 

3  3|  3  dtj  3 

SE  5?  +  S  ^ 

These  relations  need  to  be  substituted  into  equations  (2.3)  and  (2-5) 
only  for  the  derivatives  of  p  and  v  since  the  other  variables  are 
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not  associated  with  the  two  characteristic  surfaces  described  by  equation 
(4.17).  We  now  specify  that  the  curves  6  =  const,  be  tangent  to  one 
family  of  charact eristic  surfaces  and  that  the  curves  t}  =  const,  be 
tangent  to  the  other.  Hence,  the  slopes  of  the  curves  5  -  const,  and 
T)  =  const,  are  given  by  and  respectively.  Therefore,  on 

6  =  const. 

r '  -  fr =  ■  V1+QCn>  <>>-i8a) 

n 

and  on  q  =  const. 

ir ■ '  It*  -  M1+«n>  (*-i&) 

*n 

If  the  two  equations  are  applied  along  the  curve  I  =  const,  and  equation 
(4.l8a)  is  introduced,  they  can  be  combined  so  that  only  ^-derivatives 
of  p  and  v  appear: 

*tj  +  6  ^  =  1*  /P?  u  -  g  ((^)B  +  v(f)n)]gl  on.. ft  =  const. 

where  a(s,n)  =  (r  +  ensin  0)^  and  P(ti)  =  1+ e*n.  This  equation 
can  be  re-written  in  the  form 

(4.19a) 

§1  (XlP(?i,n  ‘  T(f  on  1  -  cori8t- 

which  is  the  form  given  in  reference  4o  .  Similarly,  if  the  equations 
are  applied  along  the  curve  q  =  const,  and  equation  (4.l8b)  is  used, 
they  can  be  combined  to  form 
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(4.19b) 


p5  "  £^PP  v5  =  +  ["*^pp  u  + 

[\j0(t!r)  “  v(^)  ))  77  on  q  =  const. 

Ctu  2  '  T  n  'T  n  d|  ' 

The  value  of  the  method  of  characteristics  (and  of  equations  (4.19)) 
is  that  the  equations  are  reduced  to  ordinary  differential  equations 
along  the  characteristic  curves  (in  equations  (4.19)  this  is  achieved 
only  with  the  derivatives  of  p  and  v  —  the  variables  u  and  T 
remain  as  partial  derivatives).  The  characteristic  curves  themselves 
are  functions  of  the  variables,  however,  and  must  be  determined  as  part 
of  the  solution.  The  standard  means  of  handling  this  is  by  iterative 
methods:  an  approximation  to  the  characteristics  yields  a  solution  to 
the  variables  that  can  be  used  to  obtain  a  better  description  of  the 
characteristics,  etc. 

We  now  consider  the  manner  in  which  equations  (4.19)  can  be  incor¬ 
porated  into  the  basic  finite -difference  scheme.  Consider  an  arbitrary 
mesh  point,  (m+l,N),  of  the  grid  that  was  depicted  in  figure  4.1, 
and  denote  it  with  the  subscript  3.  The  two  characteristic  curves, 

I  =  const,  and  q  =  const.,  which  pass  through  this  point  5  intersect 
the  line  s  =  mAe  at  points  that  are  denoted  by  the  subscripts  1 
and  2,  respectively  (see  figure  4.!?).  The  data  which  are  used  to 
linearize  the  finite-difference  equations  can  be  usee'  to  determine 
the  locations  of  points  1  and  2  (if  desired,  subsequent  iterations 
can  be  used  to  improve  these  values).  Points  1  and  2,  of  course,  will 
generally  not  be  grid  points,  and  evaluation  of  the  variables  there 
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will  require  the  use  of  interpolation  formulas. 


Figure  4.5-  Finite -difference  grid  for  the  characteristic  curves 


If  we  replace  the  derivatives  with  respect  to  i  and  q  with  the 
difference  quotients 


and 


const. 


dF 


const. 


and  if  we  replace  the  functions  /pp*  ,  ,  etc.  with  suitable 
averages  over  the  intervals  3-1  and  3-2,  then  equations  (4.19) 
can  be  used  to  obtain  relations  of  the  form. 


Pm+1,N  Ap 


+  B 


p  m+l,N 


D  T 
p  m+l,N 


iOa) 


(4.2 


(4.20b) 


v  —A  +  R  u  +DT 

m+l,N  v  v  m+l,N  v  m+l,N 

The  coefficients  of  these  equations  are  functions  of  the  variables  at 
points  1  and.  2  and/or  functions  of  average  values  over  the  arcs  5-1 
and  3-2.  Hence,  they  can  be  assigned  numerical  values  which  are  consis¬ 
tent  with  the  linearization  of  the  finite-difference  scheme.  Note  also 
that  the  values  of  A£  and  Aq  need  not  be  known;  only  As  remains 
in  the  equations. 

When  the  tangential  momentum  equation  and  the  energy  equation 

were  replaced  with  the  difference  equations  (4.3),  the  variable  ^ 

entered  into  the  coefficients  of  the  equations  through  the  term 

^  .  In  method  I,  an  extrapolated  value  of  p  was  used,  and  it  was 

seen  that  this  was  one  of  the  contributing  factors  to  the  occurrence 

of  instabilities.  Now,  however,  equation  (4.20a)  can  be  used  to  eliminate 

p  „  in  favor  of  u  , ,  „  and  T  ,,  this  will,  of  course,  redefine 
m+l,N  m+1,  N  m+l,N 

the  values  of  the  coefficients  of  equations  (4*3).  Once  u  and  T 
are  obtained  from  the  finite-difference  solution,  p  and  v  are 
available  from  equations  (4.20).  However,  by  making  use  of  the 
characteristic  coordinates  for  p  and  v,  the  solution  of  all  the 
variables  is  simultaneous,  whereas  in  method  I  there  was  a  sequential 
ordering  to  the  solution  cf  the  various  variables. 

Unfortunately,  for  the  blunt -body  problem  the  reduction  of  the 
differential  equations  (4.19)  to  the  algebraic  equations  (4.20)  cannot 

be  accomplished  as  easily  as  indicated  above.  Since  the  slope  of  the 

1 

u’ 

135 


characteristics  is  proportional  to 


the  slopes  are  rather  large 


near  the  axis  of  symmetry  of  the  flow  and  near  the  body  surface- 
Figure  4.6  illustrates  the  situation  by  depicting  the  two  characteristics 
which  pass  through  the  point  (s,n)  =  (0.10,  0.6o)  for  flow  around  a 
sphere  with  Y  =  1.4,  M  =  10,  Re  =  100,  b  =  0.6,  etc.  The  computation 
of  the  location  of  the  characteristics  is  based  on  the  second -truncation 
results.  Also  shown  in  figure  4.6  are  the  characteristics  which  intersect 
the  boundaries  at  s  =  0.10.  If  the  reduction  of  equations  (4.19)  as 
described  above  is  to  be  valid,  the  length  of  the  arcs  3-1  and  3-2 
must  be  small.  It  can  be  seen  that  this  will  require  an  extremely 
small  step  size  in  the  s-direction,  even  if  the  truncations  were 
accurate  enough  that  the  finite -difference  scheme  could  be  started  at 
s  =  0.10.  A  step  size  of  As  =  0.001  could  probably  be  used  but  would 
result  in  rather  large  computation  times.  If  a  more  reasonable  value 
of  As  is  to  be  used,  a  more  elaborate  means  of  integrating  equations 
(4.19)  along  the  characteristics  must  be  used.  This  would  tend  to 
destroy  the  simplicity  of  the  method  and  lead  to  rather  complicated 
computational  procedures.  Hence,  the  scheme  has  not  been  further 
investigated  although  it  is  at  least  theoretically  feasible- 

However,  the  fact  that  the  pressure  and  the  normal  velocity  can 
be  treated  by  the  method  of  characteristics  is  quite  important  in  itself, 
irrespective  of  whether  the  method  is  practical.  The  existence  of  the 
characteristic  relations  for  p  and  v  is  verification  of  the 
feasibility  of  solving  the  shock-layer  equations,  (2.2)-(2*7),  as  an 
initial-boundary-value  problem  with  the  axis  of  symmetry  for  the  initial 
station.  In  addition,  the  characteristics  may  be  useful  for  the 
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application  of  the  boundary  conditions.  In  section  C,  it  was  seen  that 
the  application  of  the  shock  conditions  led  to  complications.  The  shock 
conditions  provide  values  of  each  of  the  four  variables  at  n  =  A  as 
functions  of  the  shock  angle  0  .  This  angle  is,  of  course,  an  unknown, 
and  it  might  appear  that  the  problem  is  underdetermined.  However, 
there  is  one  characteristic  curve  that  originates  in  the  "known"  flow 
field  (i*e.,  s  <  (m+l)As  'nd  0  <  n  <  A)  and  that  intersects  the 
shock  at  s  =  (m+l)As  .  The  differential  equation  that  applies  along 
this  characteristic  provides  an  additional  constraint  on  the  values 
of  the  variables  and  permits  the  determination  of  0  .  At  the  body 
surface  there  are  boundary  conditions  on  u  (equation  (2.10)),  T 
(equation  (2.1l)),  and  v  (equation  (2.12)),  but  none  on  the  pressure 
(or  density).  The  fact  that  one  of  the  characteristic  curves  that 
originates  in  the  "known"  flow  field  intersects  the  wall  and  provides 
an  additional  constraint  indicates  that  no  boundary  condition  is  needed.'1' 
It  has,  of  course,  been  tacitly  assumed  throughout  this  investigation 
that  the  physical  boundary  conditions  were  mathematically  adequate  for 
a  solution  of  the  type  that  has  been  sought.  One  can  interprcte  the 
above  discussion  as  an  indication  that  the  boundary  conditions  are  of 
sufficient  number  to  yield  a  well  posed  initial-boundary -value  problem. 


A  discussion  of  the  use  of  characteristics  to  determine  the  number  of 
boundary  conditions  and  initial  conditions  for  partial  differential 
equations  of  various  types  is  contained  in  reference  23 • 
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In  addition,  the  use  of  the  characteristic  curves  that  intersect  the 
boundaries  may  be  a  practical  means  of  evaluating  some  of  the  boundary 
values  although  this  point  has  not  been  investigated - 

E.  Method  III 

Since  it  is  somewhat  inconvenient  to  make  explicit  use  of  the 

characteristic  equations,  we  now  consider  whether  the  basic  finite- 

difference  scheme  itself  is  capable  of  handling  the  characteristics 

implicitly.  The  use  of  explicit  finite -difference  schemes  to  solve 

simple  partial  differential  equations  of  the  hyperbolic  type  is 

discussed  in  most  texts  on  numerical  methods .  The  well-known  result 

is  that  the  domain  of  dependence  of  the  numerical  set  'me  must  contain 

the  domain  of  dependence  of  the  differential  equations.  Thus,  for 

explicit  fin1 , e-difference  methods,  a  strict  requirement  is  placed 

As 

on  the  ratio  of  the  step  sizes,  —  ,  if  the  computations  are  to  be 

stable.  However,  for  implicit  finite-difference  schemes,  no  such 
As 

requirement  on  ^  is  needed  since  the  numerical  domain  of  dependence 
extends  across  the  entire  flow  field.  This  does  not  mean  that  an 
implicit  finite-difference  method  will  necessarily  work,  but  it  does 
mean  that  the  scheme  at  least  meets  the  minimum  requirements  for 
stability.  Hence  we  now  consider  a  computational  method  in  which  all 
the  flow  variables  are  programmed  into  the  implicit  finite-difference 
scheme. 

This  method  was  briefly  examined  at  an  early  stage  of  the  present 
investigation,  and  a  rather  peculiar  type  of  instability  was  encountered 
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Although  this  method,  was  then  dropped  in  favor  of  method  I,  the  analysis 
of  the  previous  sections  has  indicated  that  this  method  may  offer  the 
hest  chance  of  success.  Hence,  the  method  is  described  in  this  section 
in  order  to  determine  the  nature  of  the  instability  that  was  previously 
encountered.  It  will  be  found  that  the  instability  is  not  related  to 
the  instabilities  encountered  in  section  C  and  that  it  can  be  eliminated 
in  a  simple  manner.  This  method  is  in  many  respects  the  simplest  and 
most  straightforward  of  the  three  methods  considered  in  this  chapter 
since  the  entire  solution  is  given  by  the  difference  relations  (4.6)- 
(4.8).  However,  there  are  certain  features  of  the  computation  that 
require  close  examination,  and  these  details  are  explained  below. 

The  partial  differential  equations  (2.2)-(2.5)  are  reduced  to 
the  difference  equations  (4.4)  by  introducing  the  difference  quotients 
given  in  equations  (4.2).  Either  the  density  or  the  pressure  may  be 
considered  to  be  the  fourth  dependent  variable  in  addition  to  u,  v 
and  T.  The  coefficients  (A,  B,  C,  and  D)  of  the  difference  equation 
(4.4)  now  contain  a  total  of  fifty-two  components. 

The  determination  of  the  boundary  conditions  (4-5)  requires  some 
explanation.  Equation  (4.5a)  contains  boundary  conditions  on  all  the 
variables,  and,  of  course,  these  constraints  are  required  by  the 
numerical  method.  Two  of  these  constraints  are  provided  by  the  slip 
and  temperature -jump  conditions,  (2.10)  and  (2.11),  and  a  third  is 
obtained  quite  simply  from  the  boundary  condition  on  the  normal  velocity, 
equation  (2.12).  Although  required  by  the  numerical  method,  there  is 
no  physical  boundary  condition  on  the  density  (or  pressure)  at  the  wall. 
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Further,  it  was  seen  in  the  preceding  section  that  the  differential 
equations  do  not  mathematically  need  such  a  constraint.  Hence,  the 
necessary  numerical  condition  must  he  obtained  from  the  differential 
equations  themselves.  It  was  noted  in  the  preceding  section  that  this 
could  be  obtained  from  the  differential  equation  that  applies  along 
the  characteristic  intersecting  the  wall.  A  simpler  means  is  to 
evaluate  the  continuity  equation  (2.2)  at  n  =  0  (note  that  the  differen¬ 
tial  equations  had  been  previously  reduced  to  difference  equations  only 

at  N  =■  1,2,... ,N  --  hence  this  condition  will  be  an  independent 
s 

equation).  The  central-difference  formulas  of  equation  (4.2)  cannot 
be  used  for  this  purpose,  but  a  forward-difference  relation  such  as 


dF 


m+l,H 


=  _  3F(0)-4f(1)+F(2)  +  0(a2) 

OAn  '  ' 


must  be  used.  The  resulting  three-point  equation  is  reduced  to  the 
required  two-point  form  as  described  in  section  C. 

A  similar  situation  arises  at  the  shock.  Although  the  Rankine- 
Hugoniot  conditions  supply  all  the  necessary  conditions  in  equation 
(4.5b),  it  does  so  in  terms  of  an  additional  unknown,  A.  Hence,  an 
additional  condition,  which  must  be  supplied  by  the  equations  themselves, 
is  needed.  The  use  of  the  characteristic  curve  that  intersects  the  shock 
was  described  in  section  D  as  one  possibility.  The  use  of  the  integral 
form  of  mass  conservation  as  described  in  section  C  is  another  possibility 
although  it  is  not  known  whether  this  will  result  in  a  numerically 
independent  condition.  Or  one  of  the  differential  equations  can  be 
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evaluated  in  such  a  way  that  a  numerically  independent  condition  is 
obtained,  as  was  done  at  the  wall. 

As  noted  earlier,  some  rough  calculations  based  on  this  method 
led  to  a  rather  strange  instability.  At  the  first  computational  step, 
the  values  of  the  variables  at  the  even-numbered  grid  points  appeared 
to  become  uncoupled  from  the  values  at  the  odd-numbered  grid  points. 
This  is  illustrated  in  figure  4.7  by  a  plot  of  the  tangential  velocity, 
u,  at  s  =  0.006.  Other  variables  show  a  similar  effect.  Unlike  the 
case  of  method  I,  this  instability  does  not  become  progressively  worse 
as  s  increases  but  is  inherent  to  the  solution  at  each  step.  It 
should  be  emphasized  that  these  results  are  accurate  solutions  of 
the  difference  equations  (4.4),  and  thus  the  source  of  this  problem 
should  be  sought  in  the  formulation  of  the  difference  equations  (and 
not  in  the  boundary  conditions  or  the  method  of  solution). 

It  has  now  been  recognized  that  the  probable  source  of  this 
difficulty  is  in  the  use  of  the  central  difference  quotient  (4.2b) 
for  the  normal  derivatives  of  p  and  v.  This  equation  does  not 
contain  the  central  point  (m+l,N),  and  this  will  obviously  have  a 
tendency  to  uncouple  the  variables  at  adjacent  grid  points.  When 
considering  the  variables  u  and  T,  this  is  no  problem  since  the 
principle  terms  involving  u  and  T  are  second  derivatives  and 
equation  (4.2c)  contains  the  variables  evaluated  at  N-l,  N,  and  N+l. 
However,  for  p  and  v,  the  first  derivatives  are  the  principle  terms, 
and  the  central  difference  of  equation  (4.2b)  is  not  appropriate. 

It  should  be  noted  that,  due  to  the  influence  of  other  terms,  the 
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coefficients  of  p(N)  and  v{N)  in  the  difference  equations  are  not 

zero;  nevertheless,  the  use  of  the  central-difference  quotient  for  the 
Sp  dv 

principle  terms,  ^  and  ^  ,  is  the  probable  source  of  the 
instability  illustrated  in  figure  4.7* 

This  conclusion  is  reinforced  by  a  consideration  of  the  order  of 
the  differential  equations  and  the  order  of  the  difference  equations. 

The  use  of  the  central  difference  quotients,  (4.2b)  and  (4.2c),  results 
in  the  difference  equations  (4.4)  being  second-order  difference  equations. 
We  have  seen  that  this  worked  well  for  the  tangential  momentum  equation 
and  the  energy  equation,  which  are  second-order  differential  equations; 
but  when  the  continuity  and  normal  momentum  equations,  which  are  first 
order,  are  included,  the  instability  of  figure  4.7  occurred.  In 
reference  Vy,  it  is  pointed  out  that  the  use  of  difference  equations 
to  approximate  differential  equations  can  lead  to  what  the  authors 
call  "computational  instabilities"  if  the  difference  equations  are 
of  higher  order  than  the  differential  equations.  These  "computational" 
instabilities  are  of  a  different  nature  than  the  more  common  instabilities 
that  usually  result  in  the  catastrophic  breakdown  of  the  computational 
procedure. 

Fortunately,  in  the  present  case,  it  would  appear  to  be  a  simple 
matter  to  remove  this  instability  by  using 

|F  .  Vl,N+l  '  +  (  s 

for  the  normal  derivatives  of  p  and  v,  although  this  point  has 
not  been  checked-  It  is  clear,  however,  that  a  successful  application 


of  this  method  must  consider  the  type  of  difference  quotients  that  are 
to  be  used. 

F.  Summary  and  Concluding  Remarks 

In  their  present  state  of  development,  the  computational  procedures 

that  have  been  presented  in  this  chapter  are  clearly  not  capable  of 

yielding  accurate  solutions  to  the  blunt-body  problem.  However,  several 

of  the  results  of  thiG  investigation  did  indicate  that  the  methods  might 

be  adequate  if  properly  formulated.  The  major  difficulties  were  traced 

to  tvo  factors:  the  determination  of  the  shock  location  and  the  handling 

of  the  terms  puv  and  pw  ,  which  appear  in  the  normal  momentum 
s  n 

equation.  In  reference  11, two  approximations  were  introduced  for  the 
constant -density  flow  in  order  to  overcome  these  difficulties.  In  the 
present  investigation  (method  i),  it  was  found  that  these  approximations 
were  less  satisfactory  for  the  variable -density  case: 

1.  The  assumption  of  a  concentric  shock  and  body  was  one  of 

the  approximations  that  permitted  the  solution  of  the  constant- 
density  flow.  The  shock  location  was  then  found  to  be  consis¬ 
tent  with  the  assumption.  However,  it  is  known  from  the  series- 
truncation  analysis  that  such  an  assumption  is  not  generally 
valid  for  a  flow  of  variable  density.  Hence,  to  specify  the 
shock  angle  in  the  present  investigation,  it  was  necessary  to 
rely  on  the  truncation  results  despite  the  fact  that  the  des¬ 
cription  of  the  shock  angle  appears  to  be  among  the  least 
reliable  of  the  results  of  the  series -truncation  analysis. 
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Near  the  axis,  the  finite- difference  computations  did  lead  to 


results  that  were  consistent  with  the  assumed  shock  angle,  but 
this  could  not  be  expected  to  be  true  over  an  extended  range 
in  s.  Furthermore,  the  fact  that  the  results  are  consistent 
with  the  assumption  does  not  in  inself  prove  that  the  assumption 
is  correct.  The  ability  to  compute  the  shock  position  and  angle 
as  a  part  of  the  finite -difference  solution  is  a  critical  factor, 
and  considerable  attention  must  be  given  thi6  point  in  any 
further  analysis  of  the  type  that  was  considered  in  this  chapter. 
It  was  noted  that  this  difficulty  can  be  circumvented  by  the  use 
of  the  thin-shock-layer  model,  but  this  must  necessarily 
influence  the  accuracy  of  the  results. 


2.  The  second  approximation  that  was  used  to  obtain  the  solution 

of  the  constant-density  equations  was  to  make  use  of  the  thin- 

shock-layer  pressure  to  evaluate  the  tangential  pressure  gradient 

do 

That  is,  the  influence  of  the  terms  puv  and  p w  was 

os  s  n 

omitted.  Since  these  two  terras  could  influence  the  flow 
variables,  other  than  r,  only  through  the  ^  term  of  the 
tangential  momentum  equation,  this  approximation  was  equivalent 


It  has  recently  been  suggested  by  Davis  [  9]  that  the  problem  should 
be  formulated  in  variables  which  are  defined  by  n  =  n/S  ,  u  =  u/u(5), 
etc.  With  these  variables,  the  shock  location  and  the  values  of  the 
variables  at  the  shock  all  become  unity.  The  unknown  quantities  such 
as  2  and  u(5),  etc.  now  appear  only  in  the  differential  equations. 
According  to  Davis,  the  difficulties  associated  with  the  computation 
of  the  shock  position,  2,  can  be  removed  in  this  manner. 


to  the  thin-shock-layer  approximation,  in  which  these  terms  are 

omitted  entirely.  The  resulting  equations  are  parabolic  in 

nature  and  hence  are  <  ^ily  treated  by  finite-difference  methods. 

For  the  variable -density  case,  the  complete  omission  of  the 

terms  puv  and  pw  leads  to  somewhat  unrealistic  values 
s  n 

of  density  and  pressure  across  the  shock  layer  for  certain  body 

shapes  (including  spheres).  Hence,  a  second  approximation, 

which  was  analogous  to  the  approximation  used  in  reference  11 

but  which  would  not  encounter  the  density  distributions  of 

thin-shock-layer  theory,  was  considered.  It  was  found  that  the 

influence  of  the  terms  puv  and  pw  through  the  density 

s  n 

was  sufficiently  large  to  cause  instabilities. 

These  two  terms  are  the  source  of  the  two  characteristics  of 
equations  (2 -2 )-(2 -5 )  that  are  of  the  hyperbolic  type.  Hence, 
a  method  was  examined  that  used  the  method  of  characteristics 
to  handle  the  computation  of  the  pressure  and  the  normal 
velocity  (method  II).  It  was  found  that  the  method  of  charac¬ 
teristics  was  somewhat  inconvenient  to  apply  due  to  large  slopes 
of  the  characteristic  curves.  Finally,  a  method  in  which  all 
the  variables  were  programmed  into  the  finite-difference  scheme 
was  considered  (method  III).  It.  was  found  that  with  this  method 
care  must  be  taken  to  choose  difference  quotients  that  are 
compatible  with  the  differential  equations.  Neither  of  these 
last  two  methods  was  investigated  in  detail.  The  last  method 
in  particular  appears  to  be  worthy  of  further  investigation. 


Despite  the  general  lack  of  success  with  these  methods,  there 
were  several  results  that  indicated  that  the  methods  might  yield  valid 
results  if  they  were  properly  formulated.  In  particular,  method  I  was 
able  to  compute  values  of  A  and  0  that  were  considerably  more  accurate 
that  the  assumptions  and  data  that  formed  the  basis  of  the  computations, 
and  the  analysis  of  method  II  indicates  that  the  formulation  of  the 
blunt-body  problem  as  an  initial -boundary-value  problem  is  justified. 

There  is,  however,  one  important  problem  concerning  this  last 
point  that  has  not  yet  been  discussed  in  detail:  it  is  not  clear  what 
effect  an  upstream  influence  has  on  the  finite-difference  procedure. 

The  series -truncation  analysis  has  clearly  shown  that  there  is  an  up¬ 
stream  influence  in  the  flow  despite  the  "parabolic-hyperbolic"  nature 
of  the  equations.  This  was  exhibited  by  the  strong  influence  of  the 
downstream  wall  temperature  on  the  flow  variables  at  the  axis.  (For 
non-spherical  bodies,  the  coefficient  Kg  of  an  expansion 
k(0)  =  1  +  KgSin  <?+•••  would  have  a  similar,  though  less  pronounced, 
effect).  Thus,  as  a  result  of  the  upstream  influence,  the  initial  data 
will  be,  to  some  extent,  incompatible  with  some  of  the  downstream 
boundary  parameters,  especially  the  wall  temperature*  We  have  seen 
that  this  incompatibility  may  be  fairly  large  if  the  first -truncation 
results  are  used  as  initial  data;  with  higher  truncations,  the  initial 
data  will  more  closely  approximate  the  data  which  are  correct  for  the 
desired  wall  temperature-  Due  to  this  incompatibility,  it  is  necessary 
to  inquire  whether  the  problem  can  be  formulated  with  a  specified 
wall  temperature.  If  this  formulation  is  to  be  valid,  the 
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finite -difference  results  must  be  able  to  quickly  adjust  to  the  imposed 
conditions  irrespective  of  the  possibly  erroneous  initial  data.  This 
is  true  of  boundary-layer  computations  —  the  influence  of  the  initial 
data  diminishes  rapidly  as  the  computations  proceed  away  from  the 
initial  station  --  but  is  less  likely  to  be  true  of  the  shock-layer 
equations  because  the  equations  have  some  hyperbolic  characteristics. 
That  is,  it  i6  likely  that  the  effect  of  the  "erroneous"  initial  data 
will  be  propagated  downstream.  Since  the  initial  data  may  not  be  valid 
for  the  particular  temperature  that  is  desired,  it  is  possible  that  a 
solution  cannot  be  found  from  the  finite-differenc  =;  procedure  unless 
the  downstream  wall  temperature  is  adjusted  to  be  compatible  with  the 
initial  data.  The  manner  in  which  the  compatible  wall  temperature 
would  be  determined,  as  well  as  the  manner  by  which  the  incompatibility 
would  be  manifested,  is  unknown  since  this  discussion  is  somewhat 
speculative. 

An  example  of  how  it  might  be  determined  can  be  seen  in  an  examina¬ 
tion  of  method  I.  It  was  there  noted  that  the  use  of  the  integral 
equation  expressing  the  overall  mass  conservation  should  automatically 
insure  that  the  Rankine-Hugoniot  condition  on  the  normal  velocity  is 
satisfied.  In  practice,  it  was  found  that  this  was  not  always  the  case. 
Normally,  the  error  was  less  than  one  percent,  but  as  instabilities 
developed,  the  error  grew,  becoming  5-10  percent  just  before  the 
computational  procedure  had  to  terminate.  This  error  possibly  could 
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have  been  removed  by  adjusting  the  value  of  the  wall  temperature  at 
each  E-step.- 

The  effectiveness  of  any  such  computation  procedure  can  be  evaluated 
by  using  the  first  truncation  results  as  initial  data.  The  wall 
temperature  that  would  be  calculated  by  the  finite -difference  scheme 
should  then  agree  to  0(sin  s)  with  the  second-truncation  temperature 
distribution  that  produced  a  concentric  shock  and  body  (Chapter  III, 
section  D.3)-  Alternatively,  if  a  constant  wall  temperature  is  specified 
together  with  the  first-truncation  initial  data,  the  results  must  make 
some  rather  obvious  and  abrupt  adjustments  (in  A  ar.d  t,  for 
example)  if  they  are  to  be  meaningful.  If  either  of  these  two  cases 
can  be  realized,  an  evaluation  of  the  reliability  of  the  r>  suits  of  a 
computation  scheme  of  the  type  that  was  considered  in  this  chapter 
will  not  be  difficult. 
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Chapter  V 


SUMMARY 

A  simplified  form  of  the  Navier-Stokes  equations  has  been  used  to 
describe  the  flow  in  the  shock  layer  of  a  blunt  body.  These  equations 
are  uniformly  valid  to  0(e )  throughout  the  entire  shock  layer  and 
thus  are  valid  in  the  same  flow  re  ime  as  the  second-order  boundary- 
layer  theory . 

The  nature  of  these  shock- layer  equations  has  led  to  the  treatment 
of  the  problem  as  an  initial-boundary-value  problem  --  with  the  axis 
of  symmetry  as  the  initial  station.  Solutions  at  the  axis  of  symmetry 
have  been  obtained  for  a  wide  range  of  flow  conditions  by  use  of  the 
method  of  series  truncation.  The  second-truncation  solutions  have  been 
found  to  yield  excellent  results  for  the  heat  transfer  at  the  wall  and 
results  for  wall  shear  that  contain  only  moderate  errors.  The  first- 
and  second-truncation  results  have  been  compared  to  the  results  of 
several  previous  investigations  of  the  flow  in  the  stagnation  region. 
These  comparisons  have  shown 

1)  that  the  basic  flow  model  is  adequate  for  values  of  the  shock 

2 

Reynolds  number  down  to  the  order  of  10  and 

2)  that,  contrary  to  the  conclusion  of  a  previous  investigation, 

the  method  of  local  similarity  may  lead  t  substantial  errors  at  the 
axis . 

With  regard  to  this  second  result,  the  ie  f  the  wall  temperature  in 
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determining  the  relative  geometry  of  the  body  and  shock  has  been  clearly 
pointed  out. 

Several  methods  based  on  the  use  of  finite -differences  have  been 
examined  as  means  of  extending  the  solutions  beyond  the  stagnation 
region.  A  method  that  was  previously  used  to  solve  the  constant- 
density  shock  layer  was  found  to  be  less  adequate  for  the  variable- 
density  case.  Although  several  interesting  results  were  obtained,  the 
method  encountered  serious  instabilities.  Two  additional  methods, 
which  should  avoid  these  instabilities,  have  been  described.  Exploratory 
computations  have  illustrated  the  features  of  each  that  require  further 
attention.  Finally,  the  results  of  this  investigation  indicated  that 
the  formulation  of  the  problem  as  an  in  tial- boundary-value  problem  is 
valid  but  that  special  attention  must  be  given  to  a  possible  incompati¬ 
bility  between  the  initial  data  and  the  vail  temperature  downstream  of 
the  axis  of  symmetry. 
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APPENDIX  A 


Coefficients  of  the  finite-difference  equations  (method  i) 

Equations  (4.3a)  and  (4.3b)  are  the  difference  equations  that  are 
used  in  method  I  (section  B,  Chapter  IV).  The  coefficients  of  these 
equations  are  determined  by  substituting  the  difference  quotients  given 
by  equations  (4.2)  into  the  differential  equations  (2.3)  and  (2-5). 

We  introduce  the  following  notation.  A  subscript  e  indicates  that  a 
quantity  is  evaluated  at  the  point  (s,n)  =  ((m+ljAs,  Mn)  but  is 
obtained  either  from  the  extrapolation  formula  (4.2e)  or  from  a  previous 
iteration.  Multiplying  the  two  equations  by  4An  and  letting  L  = 
we  obtain  the  following  relations. 


ao* -  -2[*(* -  <pv>e]  - . 


(A-l) 


B1N  1  +<€  AnN  ^pU^e  +  1 


±S**£^t  pv)  +& 
+K£AnN  v  Je  An 


(A-2) 


-A1  -  Si 

N  Ai 


,  Um,N+l  ~Um,N-l  _ 
1  An 


(A -4) 


e1n  =  0 


(A-5) 


(A-6) 
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APPENDIX  B 


Derivation  of  the  difference  equations  for  the  mass-flow  integral 
We  wish  to  show  that  the  integral  given  by  equation  *4.l4a), 

pu(r+€n  sin  0)^dn  ,  (B-l) 

can  be  evaluated  from  the  simple  algebraic  equation 

I(N)  =  A  (N)  +  Bt(n)u(n)  +  CI(N)T(N)  (B-2) 

where  the  coefficients  A^,  B^,  and  Cj.  can  be  computed  along  with 
H{N)  and  K(N)  (i.e.,  they  can  be  obtained  from  recursive  relations 
starting  at  N  =  0,  prior  to  the  calculation  of  u  and  T).  We  consider 
only  the  case  of  i  =  2  here,  but  the  more  general  case  of  i  =  4  can 
be  handled  in  an  analogous  manner. 

For  1=2,  the  difference  equation  (4.6)  may  be  written  out 
explicitly  as 

u(N)  =  H^N)  +  K  (N)u(N+l)  +  K^(N)T(N+l)  (B-3a) 

and 

T(N)  =  Hg(N)  +  K^nMn+I)  +  K^NMN+l)  (B-3b) 

where  H-^,  Hg,  K.^,  K^,  K^,  are  the  components  of  H(N)  and  K(N). 

The  integral  I  is  evaluated  by  using  the  trapezoidal  integration 
formula: 


I(N) 


2J£ 


L 
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I(N+1)-  l(N)  =  2J-1Ane  {p(N+l)u(N+l)[r+ eAn  (N+..)sin0]^ 

+  p(N)u(N)(r+€AnN  sin  0]J}  (B-4) 

Let  R(N)  =  2^  ^An  e  p(N)(r+eAn  N  sin  0)^.  If  p  is  evaluated  by 
extrapolation,  which  is  consistent  with  the  linearization  of  the  differ¬ 
ence  equations,  then  R(n)  can  immediately  be  evaluated  at  each  value 
of  N.  Equation  (b-4)  can  be  written  as 

I(N+1)  =  I(N)  +  R(N+l)u(N+l)  +  R(N)u(N)  (B-5 ) 

At  the  wall,  l(0)  =  0.  Note  that  by  recursive  application  of  equation 
(B-5)  for  N  =  1,2,...  and  by  using  equations  (B-3)  to  eliminate  the 
undesired  values  of  u  and  T,  the  integral  l(N+l)  can  be  expressed 
as  a  linear  function  of  u(N+l)  and  T(N+l),  equation  (B-2).  The 
coefficients  of  (B-2)  are  evaluated  as  follows. 

Substitute  (B-2)  into  {.  -5): 

l(N+l)  =  a^(n)  +  R(N+l)u(N+l)  +  [BT(N)  +  R(N)]u(N) 

j-  X 

+  CI(N)T(N) 

Use  equations  (B-3)  to  eliminate  u(N)  and  T(N): 

i(n+i)  =  ai(n)  +  (bi(n)  +  r(n)]h1(n)  +  ci(n)h2(n)  +{r(n+i) 

+  [Br(N)+  R(N)]Kn(N)  +  Cr(N)K21(N)}u(N+l) 

+  { [bi(n) +  r(n)]k12(n)+  ci(n)k22(n)}t(n+i) 

This  result  now  gives 
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AI(N+l)  =  A];(N)  +  [BX(N)+  R(N)J  HX(N)  +  C];(n)H2(n) 


(B-6a) 


BI(N+1)  =  R(N+1)+  [Bx(n)  +  R(n)J  ^(Nj  +  CjtNjKg^N)  (B-6b) 

and 

CjU+l)  =  [BI(N)  +  R(N)]K12(N)  +  CI(N)K22(N)  (B-6c) 

Since  Ax(o)  =  B^( 0)  =  0^(0)  =  0,  the  coefficients  A^.,  B^,  and 

can  be  evaluated  at  each  value  of  N  by  recursive  application  of 

equations  (B-6)  for  N  =  0,1,2, ...,N  • 

s 
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